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Preface 


The  work  of  this  thesis  was  prompted  by  a  paper  pub¬ 
lished  in  1975  by  J.  C.  Wyant  on  the  subject  of  phase  esti¬ 
mation  for  wavefront  correction  systems.  His  technique  for 
estimating  the  phase  differentials  measured  from  lateral 
shearing  interferometry  on  an  optical  wavefront  had  the 
form  of  a  simple  and  appealing  algorithm.  The  need  for  a 
closer  examination  of  Wyant's  technique  from  a  statistical 
point  of  view  was  a  motivating  factor  for  the  work  of  this 
thesis,  as  was  also  the  possibility  of  obtaining  improved 
performance  from  joint  processing  of  multiple  measurements 
performed  by  real-time  wavefront  correction  systems. 

I  would  like  to  thank  Capt.  Stanley  R.  Robinson  for 
suggesting  the  topic  and  giving  me  initial  advice  and 
encouragement.  I  would  also  like  to  thank  the  members  of 
my  thesis  committee,  Dr.  Peter  S.  Maybeck  and  Lt  Col 
Ronald  J.  Carpinella,  for  accepting  reading  of  this  paper 
and  for  offering  their  comments  and  constructive  criticism. 
I  also  appreciate  useful  conversations  with  2Lt  Martin  B. 
Mark  who  worked  on  a  related  topic.  I  also  thank  M.  Sally 
Lindsay  for  her  typing  of  this  thesis. 

Finally,  I  would  like  to  extend  my  gratitude  to  my 
thesis  advisor,  Capt.  James  M.  Moore,  whose  guidance  and 
encouraging  advice  throughout  this  work  are  greatly  appre¬ 
ciated. 
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Abstract 

9  Wavefront  estimation  from  shearing  interferometry 

measurements  is  considered  in  detail.  Two  analyses  are 
presented,  which  involve  the  estimation  of  constant  phase 
from  single  detector  and  detector  array  measurements.  The 
single  detector  analysis  is  carried  out  in  a  discrete  mode 
to  obtain  algorithms  based  on  photon  counting  as  the  alter¬ 
nate  means  for  use  under  low  light  level  conditions.  The 
method  used  follows  the  Maximum  A  Posteriori  and  Maximum 
Likelihood  estimation  theories.  This  is  done  for  measure¬ 
ments  made  in  both  white  Gaussian  noise  and  Poisson  shot 
noise  limited  conditions.  The  results  so  obtained  are 
trigonometric  relationships  between  the  phases  and  the 
photon  counts.  The  theoretical  performance  results  show  a 
strong  signal-to-noise  ratio  dependence.  Simulation  results 
show  that  signal-to-noise  ratios  of  17  dB  or  better  are 
needed  to  produce  adequate  estimates.  Both  theory  and  sim¬ 
ulation  show  that  an  estimate  improvement  is  obtained  as 
more  photon  counts  are  performed,  and  in  the  limiting  case, 
the  ideal  form  is  a  current  measurement.  In  this  sense, 
although  photon  counting  seems  to  be  inferior  to  current 
measuring,  the  error  variance  is  only  1.65  dB  larger  in  the 
worst  case,  where  three  photon  counts  are  performed. 

The  ML  estimator  was  found  to  be  computationally  sim¬ 
pler  than  the  MAP  estimator,  and  with  similar  performance  for 
SNR's  in  the  order  of  10  dB  and  higher. 


xiv 


An  extension  of  the  single  detector  analysis  is  made, 
using  only  the  Gaussian  noise  assumption,  to  derive  an 
algorithm  that  jointly  estimates  the  phase  distribution 
over  an  optical  wavefront.  The  procedure  is  based  on  a 
parametric  dependence  between  the  measurements  performed  by 
adjacent  detectors,  and  on  the  a  priori  knowledge  available 
through  a  covariance  matrix.  An  algorithm  for  processing 
continuous  waveform  measurements  is  developed,  but  no  com¬ 
puter  simulation  is  included  due  to  difficulties  encount¬ 
ered  in  solving  the  feedback  system  equations. 
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PHASE  ESTIMATION  TECHNIQUES  FOR  ACTIVE 


OPTICS  SYSTEMS  USED  IN  REAL-TIME 
WAVEFRONT  RECONSTRUCTION 

I  Introduction 

Reconstruction  of  a  wavefront  in  real  time  is  of  par¬ 
ticular  interest  to  the  Air  Force  because  of  the  need  to 
compensate  for  atmospheric  disturbances  and  target  varia¬ 
tions  that  adversely  affect  laser  weapons  systems.  Wave- 
front  correction  systems  of  diverse  complexity  are  employed 
to  maximize  the  irradiance  of  the  laser  on  a  target.  The 
laser  beam  is  continuously  shaped  in  real  time  by  means  of 
mirrors  to  reconstruct  the  detected  wavefront  of  the  tar¬ 
get's  radiation.  Actually,  the  complex  conjugate  field  is 
reconstructed  to  propagate  back  to  the  target  a  wavefront 
with  the  same  characteristics  but  in  complementary  form. 
Prior  to  such  reconstruction,  the  phase  distribution  of  the 
wavefront  must  be  estimated  over  the  region  of  space 
enclosed  by  the  aperture  of  the  receiving  system.  The  most 
common  method  used  for  measuring  the  phasefront  is  shearing 
interferometry.  The  search  for  improved  phase  estimation 
techniques  using  the  outputs  of  shearing  interferometers 
constitutes  the  basis  of  this  thesis.  The  shearing  inter¬ 
ferometer  will  be  discussed  in  Chapter  II. 

System  Description 

Active  optics  systems  have  been  widely  described  in 
the  literature  and  only  a  brief  description  is  necessary 
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for  the  purpose  of  this  paper.  Detailed  systems  descrip¬ 
tions  are  given  by  Hardy  (Ref  3)  ,  Hudgin  (Ref  5)  ,  Rinuner 
(Ref  10)  and  Martoni  (Ref  7:1)  among  many.  The  system 
operation  basically  consists  of  wavefront  detection,  phase- 
front  estimation,  and  beam  control.  Figure  1  shows  a  sim¬ 
plified  block  diagram  of  a  typical  system.  In  such  a  sys¬ 
tem,  a  reflecting  telescope  is  used  both  as  entrance  aper¬ 
ture  for  the  optical  radiation  from  the  target,  and  as  exit 
aperture  for  the  laser  beam.  Both  input  and  output  wave- 
fronts  travel  the  same  path  in  opposite  directions.  Part 
of  the  incoming  field  is  deflected  off  onto  a  phasefront 
sensor  usually  composed  of  two  shearing  interferometers. 

The  output  of  this  sensor  is  translated  into  control  com¬ 
mands  which  actuate  deformable  mirrors  off  which  the  laser 
beam  is  reflected  onto  the  target,  Since  this  is  done  in 
real  time,  the  atmosphere  induces  on  the  laser  wavefront 
the  reverse  distortion  effects  induced  on  the  detected 
field.  The  radiation  reaching  the  target  has,  therefore, 
been  adjusted  for  maximum  irradiance. 

Problem  Statement 

There  is  an  issue  expressed  by  the  Weapons  Laboratory 
that  when  the  target  radiance  is  low,  the  detectable  field 
is  not  strong  enough  to  perform  phasefront  estimation  base 
on  con' inuous  signal  measurements.  A  phase  estimation  tech¬ 
nique  was  proposed  by  J.  C.  Wyant  in  1975  (Ref  16:2624), 
based  on  detector  processing  of  photon  counts  observed 
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Fig  1.  A  Typical  Real-Time  Wavefront  Correction  System 
(Adapted  from  Ref  7) 

j 

during  short  time  intervals.  Although  his  assumption  is 
shot  noise  limited  detection,  his  results  are  free  from 
noise  considerations,  and  the  simplicity  of  the  resulting 
algorithm  motivates  one  to  investigate  more  deeply  into  his 
technique.  The  overall  wavefront  is  obtained  from  a  mapping 
of  independent  phase  estimates  over  the  aperture  of  the  sys¬ 
tem. 


The  purpose  of  this  paper  is  to  explore  the  photon 
counting  technique  from  a  statistical  point  of  view  and  to 
examine  in  detail  the  effects  of  receiver  and  signal  shot 
noises  in  order  to  determine  the  extent  to  which  this  pro¬ 
cedure  can  be  applied  efficiently.  Joint  processing  of 
multiple  detectors  has  also  been  considered  within  the 
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scope  of  this  work  in  an  attempt  to  obtain  improved  results 
over  linear  mappings  currently  used  tc  estimate  the  wave- 
front  over  the  region  of  interest. 

Approach 

The  basic  approach  to  the  problem  is  based  on  countable 
observables  obtained  by  decomposing  the  continuous  output 
of  a  detector  into  discrete  components.  The  problem 
reduces  to  classical  parameter  estimation  theory  and  will 
be  carried  out  using  the  concepts  of  Maximum  A  Posteriori 
and  Maximum  Likelihood  estimation  theories.  To  fulfill  the 
purpose  of  this  approach,  it  is  assumed  that  specially 
designed  detectors  are  available,  based  on  the  promising 
future  of  charge-coupled  devices  (CCD)  (Ref  13 -.chapter  12), 
which  are  capable  of  integrating  the  detected  field  signal 
over  short  periods  of  time  and  dumping  the  contents  into 
the  registers  of  a  computer.  Such  a  receiver  can  be  repre¬ 
sented  mathematically  with  a  time  correlator  to  be  described 
later  in  Chapter  III. 

Scope  and  Assumptions 

In  this  thesis,  the  wavefront  process  will  be  con¬ 
sidered  slow  varying  in  time  such  that  a  stepwise  approxi¬ 
mation  to  the  actual  variation  can  be  performed.  The  phase 
in  each  step  of  length  T  will  be  constant  and  the  analysis 
of  the  problem  will  be  limited  to  a  single  observation 
interval  (0,T) .  An  extension  to  sequential  estimation  over 
successive  intervals  can  then  be  performed  using  Gauss- 
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Markov  parameter  models.  That  will  not  be  considered  in 
this  thesis.  Constant  phase  is,  therefore,  the  basic 
assumption  of  the  forthcoming  developments.  The  disturb¬ 
ance  induced  by  the  atmosphere  will  be  modeled  as  an  addi¬ 
tive  noise  phase  to  the  target  wavefront.  There  is  no  need 
to  distinguish  between  the  target  and  noise  induced  phases 
because  the  reverse  disturbance  effects  produced  on  the 
laser  wavefront  by  the  atmosphere  cancel  the  noise  compon¬ 
ent.  Therefore,  they  will  be  lumped  together  into  a  single 
parameter  0,  where  0  is  a  random  variable.  The  probabilis¬ 
tic  descriptions  of  0  will  be  fitted  to  the  ones  given  by 
Gaussian  and  uniform  probability  density  functions  in  the 
interval  (— tt  ,  tt  )  .  Other  than  for  ease  in  estimator  deriva¬ 
tion,  the  Gaussian  model  is  chosen  considering  that  for 
slow  varying  wavefronts,  the  phase  variations  are  more 
likely  to  be  concentrated  about  the  zero  value  and  less 
likely  as  the  phase  value  increases.  On  the  other  hand, 
because  the  sensor  output  is  a  sinusoidal  variation,  the 
uniform  density  is  also  a  logical  choice  since  the  phases 
are  equally  likely  in  the  interval  ( — tt , tt )  . 

Estimation  of  0  will  be  analyzed  in  the  presence  of 
noise  from  two  points  of  view:  predominant  detector  noise 
and  predominant  signal  shot  noise  processes.  Chapter  IV  is 
devoted  to  the  analysis  of  the  detector  limited  case,  where 
the  noise  is  modeled  as  a  continuous  white  Gaussian  random 
process.  Chapter  V  is  devoted  to  the  analysis  in  signal 
induced  noise,  where  the  noise  is  modeled  as  a  discrete 
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Poisson  count  process.  A  mixed  mode  of  continuous  and  dis¬ 
crete  processes  will  not  be  considered  in  this  paper,  but 
deserves  future  attention.  An  extension  of  the  Gaussian 
noise  analysis  is  performed  in  Chapter  VI  where  the  same 
concepts  are  applied  to  joint  processing  of  two  plane 
detector  arrays.  Finally,  Chapter  VII  makes  a  summary  of 
results  and  conclusions,  and  presents  suggestions  for  fur¬ 
ther  study. 
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II  The  Shearing  Interferometer 


The  wavefront  sensor  of  interest  in  this  thesis  con¬ 
sists  of  two  ac  heterodyne,  lateral  shearing  interferometers. 
This  sensor  configuration  is  depicted  in  Figure  2.  The 
field  entering  the  system  is  beamsplit  into  two  channels. 

Each  channel  has  a  shearing  interferometer  composed  of  two 
confocal  lenses  which  constitute  a  Fourier  transform  pair. 

The  field  at  the  common  focal  point  is  the  Fourier  trans¬ 
form  of  the  received  field.  The  Fourier  transform  field  is 
sampled  with  a  radial  grating  displaced  off  the  optic  axis 
and  rotating  with  velocity  v  and  period  t .  An  expanded 
view  is  shown  in  Figure  3. 

If  the  aperture  is  located  at  the  front  focal  plane  of 
lense  ,  the  field  observed  at  the  back  focal  plane  of 
lense  L2  is  sheared  into  a  number  of  components  laterally 
displaced  from  each  other  by  an  equal  distance  Ms^,  known 
as  the  shear  distance,  where  s^  is  a  parameter  of  the  inter¬ 
ferometer  determined  by  v,  and  M  is  a  magnification  factor 
determined  by  the  ratio  of  the  lenses  in  the  system.  The 
field  components  so  displaced  interfere  with  each  other,  and 
the  modulated  output  is  observed  in  the  form  of  an  inter- 
ferogram  spread  over  a  detector  array.  It  is  this  slowly 
but  continuously  changing  interferogram  that  contains  the 
wavefront  phase  information  being  sought. 
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Rotating  Grating 


Fig  2.  A  Two  Shearing-Interferometer  Sensor  Used  with 

Wavefront  Correction  Systems  (Adapted  from  Ref  3) 

Output  Field 

The  sheared  field  at  the  focal  plane  takes  on  slightly 
different  forms  for  broadband  (white  light)  and  monochro¬ 
matic  fields.  However,  the  same  equation  is  applicable  in 
both  cases  when  the  shear  is  small.  Thus,  for  spatially 
coherent,  white  light  aperture  fields  when  the  shear  dis¬ 
tance  Ms^  is  small,  the  detector  field  intensity  is  given  by 
(Ref  6:58,61) 

Vtl.t1  -  (I  +  iK'A<Fa>)2 

+  —  M2[A(r  )]2cosC(?  )  sin  (ut  +  0  (r  )) 

7T  cl  Si  f  u  S  >  t 

for  0  <  t  <  T,  (2-1) 
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Fig  3.  Expanded  View  of  the  X-Channel  Interferometer 
(Adapted  from  Ref  6) 


are  phase  functions  in  terms  of  the  phases  at  locations  ±Ms^ 
from  ra»  The  term  A(ra)  is  the  amplitude  of  the  aperture 
field  at  location  r^  =  Mr^  shown  in  Figure  3.  For  slow 
varying  fields  over  the  aperture  and  constant  intensity, 

Eq.  (2-1)  can  be  simplified  to  (Ref  7:10,6:59) 
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The  frequency 

w  (2-5) 

is  the  fundamental  modulation  frequency  of  the  detector 
field  (Ref  6:41).  Double  frequency  and  higher  order  fre¬ 
quency  terms  have  been  dropped  from  Eqs.  (2-1)  and  (2-4), 
anticipating  subsequent  signal  processor  filtering. 

Detector  Signal 

The  output  r(t)  of  a  detector  at  location  r  in  the 

dl 

back  focal  plane  is  computed  from  the  received  field  intensi¬ 
ty  ^d^d  t,f^ra  t^  and  *-he  detector  noise  as  follows: 

r(t)  =  s(t,0(t))  +  n(t),  at  r^  t  Aperture,  (2-6) 

where  (Ref  2:54-55) 

s (t f 6 (t) )  -  /  ia(?a<t.e(MFd(tna?a 

O  Ad 

for  0 (r^  t)  constant  over  A^,  (2-7) 

is  the  signal  current.  The  constant  n/hf  is  a  detector  con- 

o 

version  factor,  A^  is  the  area  of  the  detector,  and 
Id^rd  t,6^rd  t^  is  the  of  the  detector  given  by  Eq. 

(2-4)  at  location  r,  =  r  /M.  The  assumption  made  here  is 

1  al 

that  the  field  intensity  at  r^  is  constant  over  the  detec¬ 
tor  area.  If  the  observation  time  is  short  enough  such  that 
the  phase  is  approximately  constant  in  the  interval  (0,T)  to 
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fit  the  assumptions  of  Chapter  I,  the  time  dependence  of 
the  phase  can  be  dropped.  The  interval  T  is  assumed  to  be 
much  larger  than  the  period  x  of  the  modulated  field. 
Equation  (2-6)  becomes 

r (t)  =  s  ( t , 0 )  +  n (t) ,  0  <  t  £  T  (2-8) 


where 


s(t,8)  =  a  +  b  sin(ot  +  0)  (2-9) 

is  the  output  signal  of  one  detector,  and 

*  -  S1-  (t  +  4 )m2a\  <2-l0> 

o  \  TT  / 

b  a  .  (2-u) 

Chapter  III  further  models  the  detector  signal  in  order  to 
pursue  the  photon  count  approach  established  in  Chapter  I, 
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Ill  Phase  Estimation  Preliminaries 


Reconstruction  of  the  wavefront  requires  knowledge  of 
the  phase  distribution  of  the  optical  field  within  the  aper¬ 
ture.  However,  the  knowledge  provided  by  the  shearing 
interferometers  is  in  the  form  of  a  phase  difference  between 
two  points  along  the  line  of  shear  as  given  by  Eq.  (2-3) . 

In  the  simplest  form,  the  problem  is  of  phase  difference 
estimation  from  the  measurement  of  a  single  detector.  The 
phase  difference  is  known  as  the  wavefront  difference  func¬ 
tion  given  by  (Refs  6:49  and  7:16) 

Ac()(ra)  =  <f>  (?a  -  Mij)  -  <{>  (ra  +  Msd)  (3-1) 

where  <J>  (r  -  Msd)  and  $  (?a  +  Msd)  are  the  phases  at  the 
aperture  which  we  ultimately  want  to  estimate.  The  compo¬ 
nents  of  the  wavefront  difference  function  in  cartesian 
coordinates  are  given  by 


<Mra  -  Msd) 

■  +  (Xa  -  Hsd*'V' 

X-shear 

(3-2) 

-  Ya  -  MEdy>' 

Y-shear 

(3-3) 

<j>  (ra  +  Mr  ) 

‘  *<Xa  +  Msdx'V' 

X-shear 

(3-4) 

-  Ya  +  Msdy>’ 

Y-shear 

(3-5) 

Figure  4  illustrates  the  phase  distribution  described  by 
Eqs.  ( 3—2 )  —  (3  —  5) .  The  shear  distances  are  usually  made 
equal  so  that  sdx  =  s^  —  s^.  When  M  =  -1  (rd  and  sd  are 
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Fig  4.  Phase  Distribution  at  the  Aperture  Seen  by 
Detectors  in  the  Back  Focal  Plane 

also  negative) ,  the  wavefront  difference  function  measured 
at  the  detector  plane  gives  information  of  the  phase  points 
located  one  shear  distance  away  on  each  side  of  the 
detector.  To  estimate  the  actual  phases  <£,  the  information 
must  be  collected  from  all  detectors.  This  will  be  addressed 
in  Chapter  VI.  Chapters  IV  and  V  will  concentrate  on  the 
estimation  of  the  phase  difference  0  (r^)  from  one  detector 
only.  Prior  to  considering  the  estimation  problem,  a  signal 
model  must  be  found  using  the  hypothetical  CCD  photon  count¬ 
ing  detector. 


< 
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Signal  Model 

A  signal  model  is  developed  in  this  section  to  estimate 
the  phase  difference  6  measured  by  each  detector  following 
the  integration  approach  with  CCD  devices.  This  is  equiva¬ 
lent  to  saying  that  the  output  of  each  detector  is  sampled 
by  a  correlator  structure  as  shown  in  Figure  5.  In  this 
correlator  model,  the  observation  r(t),  which  represents 
the  cathode  current,  is  sampled  by  a  vector  of  K  orthogonal 
basis  functions 


i(t)  =  [^(t)  ...  i|>.(t)  ...  i^k(t)]T,  (3-6) 


(the  superscript  T  means  transpose  of  the  matrix)  where  each 
basis  function  ip ^  (t)  is  given  by 


(t) 


-  ,  t.  <  t  <  t 
q  J  3+1 


=  0  ,  otherwise, 


(3-7) 


with  q  being  the  electron  charge  and  t^  a  sequence  of  k 
equal-length,  non-overlapping  time  subintervals  in  one 
observation  interval  of  length  T.  The  resulting  output  is 
a  vector  of  k  observations 

r  =  [r1  ...  r j  ...  rk)T  (3-8) 

which  represent  photon  counts  in  the  time  subinterval 
(tj,tj+l>.  the  signal  represented  in  this  manner,  the 

following  equalities  are  true: 
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r  . 
3 


Fig  5.  Equivalent  Correlator  Model  of  the  CCD  Detector 


ijj^  (t)  j  (t)  dt 


2  '  1  3 

q  k 


0  ,  i  /  j 


(3-9) 


and 


2,  k 

r  (t)  =  l.i.m.  2_r-  z  r  ,i|> .  (t)  , 

Jc  ->  «,  1  j=i  3  3 


(3-10) 


is  satisfied.  This  condition  known  as  the  Cauchy  criterion 
for  convergence  of  random  sequences  (Ref  8:262)  is  satis¬ 
fied  by  the  choice  of  the  t^'s  as  non-overlapping  time 
subintervals.  Equations  (3-9)  and  (3-10)  are  the  tools 
needed  in  the  forthcoming  parameter  estimation  analysis. 
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This  analysis  will  be  the  classical  discrete  observation 
problem  with  generalizations  to  the  continuous  observation 
form  using  Eq.  (3-10) .  For  this  purpose  and  to  define  t^ 
properly,  the  observation  period  is  divided  into  k  equal 
subintervals  (tj,t_.+1).  The  functions  ^  (t)  are,  therefore, 
identical  except  for  their  sequential  position  in  time.  In 
order  to  simplify  the  derivation  let  the  initial  time  be 
t^  =  -T/2k  and  the  final  time  t^+T  be  t^  =  (2k-l)T/2k  for 
any  interval  T.  Thus,  the  integration  limits  for  the  jth 
subinterval  are 


_  (2 j-3) T 

j  _  2k 


(3-12) 


and 


t . 


j+1 


-  ( 2  j  - 1 )  T 

2k 


The  correlation  operation 


r  (t)  \p  j  (t)  dt  , 


indicated  in  Figure  5  is  then  given  by 

r(t) 


-fcj+l 


dt 


(3-13) 


(3-14) 


f" j+1  f  3+1 

=  /  51^01  dt  +  j  at  , 


-Vi 


3  3 

which  can  be  conveniently  defined  as 


(3-15) 


rj  -  sj(8)  +  nj  ' 


(3-16) 
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where 


ftj+1 
J  sit. 


s  j  (0 )  = 


0 )  dt 


(3-17) 


n .  = 

3 


rVi 


n (t) dt 


(3-18) 


with  tj  and  t^+1  given  by  Eqs.  (3-12)  and  (3-13).  Because 

the  observation  r(t)  has  the  form  of  an  electrical  current, 

the  correlation  functions  of  Eq.  (3-7)  were  selected  so  that 

r j  ,  s  ( 0 ) ,  and  n^  have  the  physical  interpretation  of  photons, 

both  signal  and  noise,  counted  in  each  observation  subinter- 

k 

val  (tj,tj+1).  The  sum  £  r^  corresponds  to  the  total 
photons  counted  in  the  observation  interval  (t^t-^+T).  If 
the  functional  form  of  s(t,0)  given  by  Eq.  (2-9)  is  substi¬ 
tuted  into  Eq.  (3-17) ,  the  signal  in  discrete  form  becomes 


s j  (0 )  = 


(a  +  b  sin (wt  +  0)]dt  ,  (3-19) 


for  0  (t)  constant  over  T. 

The  integration  indicated  by  Eq.  (3-19)  yields 


V61  “  I  >Vi  -  V 


+  -  [cos (wt j  +  0)  -  cos(wtj+1  +0)]  .  (3-20) 


If  the  observation  interval  is  several  times  the  modulation 
period  so  that  many  frequency  cycles  are  observed  each  time, 
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then  they  are  related  by 


N  = 


T 

T 


(3-21) 


where  t  is  as  in  Eq. 
substituted  into  Eq. 


s.  (0) 


(2-5)  . 
(3-20)  , 

bT 

2rqN 


With  the  identity  of  Eq.  (3-21) 
the  signal  equation  becomes 


[cost—  (2j  — 3)+Q)— cos  (™(2j-l)+0)  ]  ,  (3-22) 


where  the  definitions  of  Eqs .  (3-12)  and  (3-13)  were  used 

as  integration  limits.  However,  in  order  to  preserve  the 
phase  information,  a  small  shear  distance  is  required  (Refs 
6:41  and  16:2622).  For  small  shear,  the  modulation  fre¬ 
quency  a)  must  be  low  relative  to  the  optical  field  frequency 
such  as  25  kilohertz  (Ref  4:363).  Thus,  typically  x  =  0.04 
milliseconds  and  N  can  range  from  unity  to  a  very  large  num¬ 
ber  depending  on  how  fast  the  phase  changes  in  time.  For 
simplicity  and  to  check  with  Wyant's  results,  let  N  =  1  so 
that  Eq.  (3-22)  becomes 


sj  (0) 


aT  bT 
q’k  2nq 


[cos(£(2j-3)+Q)-cos (£(2j-l)+G) ]  .  (3-23) 

Although  Eq.  (3-23)  gives  an  exact  expression  for  the  photon 
counts  in  terms  of  0,  it  is  not  in  a  workable  form.  It  is 
therefore  necessary  to  expand  the  cosine  functions  to  obtain 
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Sj(0)  =  qf  +  rSq"  |  [cos£(2j-3)  -cos£(2j-l)  ]  cosG 

-  [sin^-(2  j-3)  -sin^-(2  j-1)  ]  sinoj  . 

(3-24) 

Equation  (3-24)  is  the  expression  that  will  be  used  to  rep¬ 
resent  the  signal-generated  photon  counts.  For  ease  in  fur¬ 
ther  derivations,  Eq.  (3-24)  can  be  more  compactly  written 
as 

sj(e)  =  +  ^  [<*jC°se  -  8 j sinG  ]  ,  (3-25) 

where 

dj  =  cos  (2 j-3)  -  cos  ^  (2 j-1)  (3-26) 

8j  =  sin  |  (2 j-3)  -  sin  £  ( 2 j-1 )  .  (3-27) 

With  the  above  signal  model  completed,  the  estimation 

A 

procedure  to  obtain  0 (r)  will  be  considered  next.  The  "hat" 

over  0  (r)  indicates  that  this  is  just  an  estimate  of  0,  and 

the  argument  "r"  indicates  that  the  estimate  is  in  terms  of 
the  observation  vector  r  of  Eq.  (3-8).  The  estimate  6  (r) 
will  be  performed  based  on  the  criteria  of  Maximum  A  Pos¬ 
teriori  and  Maximum  Likelihood  estimation.  Before  applying 
these  criteria  to  find  the  phase  estimates  in  Chapters  IV, 

V,  and  VI,  these  concepts  will  be  briefly  explained  in  the 
following  section. 
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Maximum  A  Posteriori  and  Maximum  Likelihood  Estimation 
Theories 

Estimation  of  the  phase  difference  0  measured  by  a 
single  detector  at  a  fixed  location  r^  and  time  t,  and  of 
the  actual  phases  $  measured  over  the  aperture  by  two  detec¬ 
tor  arrays  will  be  made  applying  the  concepts  of  Maximum  A 
Posteriori  (MAP)  and  Maximum  Likelihood  (ML)  estimation 
theories.  The  MAP  and  ML  estimates  of  Q  are  those  values 
0 (r)  for  which  the  probability  of  having  found  the  correct  6 
after  the  measurement  is  made  is  maximum.  This  is  equiva¬ 
lent  to  maximizing  the  a  posteriori  probability  density 
function  of  0  conditioned  on  the  observations.  Let  this 
a  posteriori  density  be  represented  by 


f0|r(O|£)  ,  (3-28) 

and  let  it  be  maximized  by  the  proper  choice  of  0 (r) .  With 
the  condition  that  the  maximum  occurs  within  the  range  of  6 , 
maximization  is  obtained  by  setting 

^[£e|£(0|r)J  -  0  .  (3-29) 

An  equivalent  and  sometimes  more  convenient  form  of  Eq. 
(3-29)  is  given  by 


^Unf0|r<MrJ]  ■  0  . 


(3-30) 


By  using  the  Baye ' s  rule  substitution 


r(0|r)  = 


f  r  |  0  (-rl°>f0  (0> 

*7^ 


(3-31) 
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In  Eq.  (3-30),  the  necessary  but  not  sufficient  condition 
for  the  MAP  estimate  is  found  to  be  (Ref  14:58) 

^[lnfr|Q  (r|6)  +  Inf  Q  (9 )  ]  |  =  0  .  (3-32) 

MAP 

Because  0  is  in  the  argument  of  a  sine  function,  it  is 
modulo  2tt.  The  ML  estimate  of  the  phase  is  then  found  by 
modeling  the  a  priori  density  fg  (0)  as  uniform  within  the 
range  of  0,  which  is  the  interval  (-n,u).  The  ML  estimate 
is,  therefore,  given  by  (Ref  14:65) 

^[lnfr|0  (r  j  0 )  ]  |  =  0  .  (3-33) 

ML 

Equations  (3-32)  and  (3-33)  are  the  equations  for  MAP  and 
ML  parameter  estimation.  The  parameter  so  found  is  the  one 
with  highest  probability  of  being  the  true  value.  There,  of 
course,  may  be  false  solutions,  and  an  error  is  associated 
with  each  estimate.  The  errors  will  be  treated  in  the  next 
chapters  after  the  solution  algorithms  have  been  found. 

The  problem  of  estimating  0 (r  )  in  the  presence  of  noise 

Si  f  t 

will  be  addressed  next;  first,  it  will  be  investigated  in 
the  context  of  a  Gaussian  problem  where  the  only  noise  is 
due  to  thermal  noise  limited  detectors,  and  then  in  the  con¬ 
text  of  a  Poisson  problem  when  the  shot  noise  is  the  predom¬ 
inant  noise  source.  The  approach,  however,  will  be  the  same 
in  both  Gaussian  and  Poisson  cases:  photon  counting  as  the 
only  alternative  under  low  light  level  conditions. 
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IV  Phase  Estimation  in  White  Gaussian  Noise 

When  the  output  field  of  the  interferometer  is  received 
by  a  detector  array  with  predominant  thermal  noise,  the 
noise  is  adequately  modeled  by  a  zero-mean,  stationary 
white  Gaussian  process  described  by  the  probability  density 
(Ref  8:360) 

fn(t)(n)  =  ^tjo^]~1/2  exp[-n2/2a^]  ,  (4-1) 


where  the  variance  is  given  by 


(4-2) 


The  process  has  a  double-sided,  flat  power  spectral  density 
N 

over  and  beyond  the  spectral  region  of  interest.  If  the 
noise  is  described  in  the  form  of  a  random  current,  addi¬ 
tive  to  the  current  output  of  the  detector,  the  spectral 
density  is  given  by  (Ref  8:361) 

N  2kT 

=  — — -  (watts  per  hertz  —  per  ohm)  ,  (4-3) 

Z  R 

e 

where  k  is  the  Boltzmann  constant,  Tq  is  the  operating  tem¬ 
perature  and  R  is  the  equivalent  resistance  of  the  detec- 
S  No 

tor.  On  a  per-ohm  basis,  has  the  units  of  energy.  The 
correlation  function  of  the  noise  process  is  the  inverse 
Fourier  transform  of  the  power  spectral  density.  For  sta¬ 
tionary  white  noise  it  is  given  by 

N 

E [n (t) n ( t ' ) 3  =  ~y  S(t-t')  .  (4-4) 
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On  n  per-ohm  basis,  the  correlation  function  has  the  units 
of  power.  The  first  order  statistics  are  .jivon  by 

E  [ n  ( t. )  ]  =  0  .  (4-5) 


Before  applying  the  cr..cepts  of  MAP  and  ML  estimation, 
it  is  necessary  to  obtain  the  density  function  of  the  obser¬ 
vation  vector  r  conditioned  on  the  parameter  0,  as  required 
by  Eqs.  (3-32)  and  (3-33).  Therefore,  the  following  devel¬ 
opment  is  made: 

Referring  to  Eqs.  (3-14)  and  (4-6) 


Referring  to  Eqs.  (3-9)  and  (4-4) 


E[n?] 


dt'E[n(t)n(t')Hj  (t)  (t '  ) 


N 

dt'  [-£  5(t-t,)^j(t)’i;j(t')] 


•Jm (t)dt 


N  T 
o 

2q2k 


(4-7) 
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There!  re, 


E  [r  j  ]  =  Sj{0)  (4-8) 

NT 

Var  [r  .  ]  =  Var[n.]  =  — .  (4-9) 

3  3  2q  k 

In  addition  to  the  variance  and  the  mean  of  r ^ ,  it  is  neces¬ 
sary  to  demonstrate  that  r ^  and  r^  are  uncorrelated.  Thus, 

cov  [r  j  ,  r  ]  =  E  [  (r  j  -  E[r..])(rj  -  E[rg])] 

=  E[  (r  -  s.  (0)  )  (r  -  s  (0)  )  ] 

J  J  Si  4. 


Further  evaluating  Eq.  (4-10) 


(4-10) 


E‘„3nq] 


dt'EInltlnlt')  J*.  (t>*  Ct  ’  >  J 

J  SI 


dt'  [-£  6  (t-t ' )  <j>  •  (t)  tp  (t ' )  ] 
*•  j  q 


T 

Therefore,  from  Eqs. 


cov [r . , r  ] 

1  q 


(t)^q(t)dt  =  0,  j  f  q  . 

(3-16) ,  (4-6) ,  and  (4-11) 

=  E  [n .  n  ]  =  0  ,  j  /  q  , 

J  M 


(4-11) 


(4-12) 
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Sj  (0)  ) 


(4-14) 


The  MAP  estimate  of  a  parameter  0  observed  in  noise  is 
based,  as  indicated  by  Eq.  (3-32),  on  the  conditional  dens¬ 
ity  of  the  observations  r  given  0,  as  well  as  on  the 
a  priori  density  of  0.  Considering  first  the  conditional 
density  of  r  given  0,  the  following  procedure  is  developed: 


From  Eq.  (4-14)  , 


fr|e  (-l6)  = 


'  2.  1 

aJs 

v  N 

o-l 


k/2 


exp 


2  k 

H  .E,  -  “*<6)  *  s2(0>) 

o  3  =  1  J  J  J 


(4-15) 


The  natural  logarithm  of  Eq.  (4-15)  is 


lnfr|e(-^0)  2  ln  ttN  T  N~T  ,Z  (rj  ~  2rjsj(0)  +  sj(0)) 

—  o  o  3  =  1  J  J 


=  2a.  K  £  r  s  (q  ) 
NQT  rjSi(G) 


-  £  s2  ( 0 ) 

NT  .L,  Sj 
o  3=1  J 


_  a!k  £  r2  +  —  In  -a?>L 
NT  .  -  rj  2  in  ttN  T  * 
o  3=1  J  O 

The  derivative  of  Eq.  (4-16)  with  respect  to  0  is 


(4-16) 


30  ln  fr  |  0  0) 


^  £  r  JL  s  ( Q ) 

N0T  rj  30  Sj(9) 


-  T&  Ve)  w  V61  • 

o  3=1  J  J 

(4-17) 

By  substituting  Eq.  (4-17)  into  Eq.  (3-32) ,  a  general 
expression  for  the  MAP  estimate  is  obtained  and  is  given  by 


2  k 

ifr  2  trj  "  sj(6nal  sj(6)  +  Wllnf0(O)1  *  o  • 

o  3=1 

(4-18) 

The  solution  of  Eq.  (4-18)  when  the  appropriate  expressions 
for  Sj(0)  and  fg  (0)  are  substituted  for  yields  Ofr)]^^ 


Before  doing  this,  it  is  interesting  to  observe  what  happens 
when  k  becomes  very  large.  It  is  very  simple  to  prove  that 


g(t)  =  l.i.m.  A  E  g.<Ji.(t) 
k  -*•  »  j=l  11 


=  l.i.m.  A  E  h.<}>.(t)  , 

k  -*■  00  j=l  1  1 


(4-19) 


where  (t)  are  orthogonal  functions  such  that 


/ 


it)  -  j 


i  i  ( t )  4> 


(t)  =  0  ,  i  f  j  , 


(4-20) 


then  l.i.m.  A 

k  -*■  " 


k  fT 

3-1  1  1  l 


g  (t ) h (t) dt  .  (4-21) 


With  the  help  of  Eqs.  (4-19) - (4-21)  and  the  orthogonal  rela¬ 
tion  of  Eq.  (3-10)  ,  Eq.  (4-18)  becomes  (t^-*0  as  k-*00) 


~  I  (r  (t)  -  s(t,0)]^-  s(t,0)dt  +  ^[lnfe(8)]  =  0  . 

o  J~ 


(4-22) 


Equation  (4-22)  is  the  general  expression  for  the  MAP  esti¬ 
mate  0  of  a  parameter  0  in  a  continuous  waveform  r(t) 
observed  in  white  Gaussian  noise  (Ref  14:275). 

Returning  to  Eq.  (4-18),  the  substitution  for  Sj(0) 

A 

given  by  Eq.  (3-25)  is  made  to  obtain  0 (r)  as  follows: 
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From  Eq.  (3-25) 


;a  rp  Kryi 

s.(0)  =  —r  +  ■»— - ■  [ci.cosO  -  g.sinO] 

3  qk  2?rq  3  3 


(4-23) 


Its  phase  derivative  is 


3  e  ,n.  _  _bT 

30  sj(G)  2TTq  l“j 


[a^sin0  +  BjCOs0] 


(4-24) 


Equation  (4-18)  can  now  be  written  as 


T?TT  J1  J(rj  -  f£  -  25}  l“jcose  -  6jsin0l> 

-  (a j sin©  +  $jCos0])j 


+  ^  [Inf Q  (0 ) ]  =  0  . 


(4-25) 


In  order  to  simplify  notation,  let  a  new  function  F(0)  be 
defined  as 


F ( 0 )  -  jq  [lnf9  (0))  . 


(4-26) 


After  some  algebraic  manipulation  of  Eq.  (4-25) ,  the  follow¬ 
ing  MAP  estimate  of  0  is  obtained: 


aT, 


I  |[r.  -  ^)a.sin§(r)  +  [r .  -  ^J6.cos0(r); 
j=1|  3  qk  3  -  3  qkJ  3  - 1 


bx 


2-q  jijIV  2 


)si 


^-Jsin20  (r)  +  a.0.cos20(r) 


3  3 


ttN 
_ o 

qkb 


F  (0 ) 


0  , 


(4-27) 
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whore  u_.  and  3_.  arc  defined  by  Eqs.  (3-26)  and  (3-27). 
Although  Eq.  (4-27)  defines  as  is,  it  is  not  in 

final  form.  The  second  summation  in  Eq.  (4-27)  vanishes 
for  all  values  of  k  except  for  k  =  2.  This  was  determined 
on  a  computer  check  and  will  be  used  without  a  rigorous 
proof.  Therefore,  Eq.  (4-27)  becomes 


k 

£ 

j  =  l 


— ]  a  .  sinO  ( r )  +  £ 

qk  3  -  3=1 


1^1 6 jCOse (r) 


ttN 

"  ^  F<0>  =  0  ;  k  >  3  ,  (4-28) 


whore  the  constraint  that  k  >  3  is  imposed  on  the  estimator 
because  for  k  =  1,  the  two  sums  vanish. 

Now,  the  conditions  for  which  the  second  summation  in 
Eq.  (4-27)  vanishes  are  to  be  investigated.  For  that  to 
happen,  it  is  required  that 


k  (a? 
^ 

j=l 


2 


0 


(4-29) 


k 

and  £  a.B-  =  0  .  (4-30) 

j  =  l  3  3 

Substitution  of  the  expressions  given  by  Eqs.  (3-26)  and 
(3-27)  for  ou  and  Bj  into  Eqs.  (4-29)  and  (4-30)  yields 

^  TT  1  2 

£  {  [cos  ?-(2  j-3 )  -  cos  ~{2j-l)]z 

j  =  l  K  K 


-[sin  £(2j-3)  -  sin  £(2j-l)]2}  =  0  ,  (4-31) 
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and 


l  {[cos  ?- ( 2 j  —  3 )  -  cos  £-(2j-l)] 
j=l  K  K 

[sin  £-(2j-3)  -  sin  ^(2j-l)]}  =  0  .  (4-32) 

By  carrying  out  the  operations  indicated  by  Eqs.  (4-31) 
(4-32) ,  it  is  found  that  the  conditions  for  a  vanishing 
term  are 

£  {cos  ^r(2j-3)  +  cos 
j=l  K  * 

-  2cos  ^(2j-2)  }  =  0  (4-33) 

and  £  {sin  ^-(2j-3)  +  sin  ^(2j-l) 

j=l  K  k 

-  2sin  ^-(2j-2)}  =  0  .  (4-34) 

The  question  is,  for  what  values  of  k  do  Eqs.  (4-33)  and 
(4-34)  hold?  The  easiest  way  to  find  out  is  by  computing 
the  summations  for  a  number  of  k's.  This  was  done  on  the 
computer  for  k  =  1  to  k  =  28  with  the  assuring  result  that 
only  for  k  =  2  both  sums  do  not  equal  zero.  It  was  also 
observed  (and  will  be  used  without  a  rigorous  proof)  that 

k  k 

£  a  =  £  6.  =  0  (4-35) 

j=l  1  j=l  D 

for  all  k. 

Therefore,  Eq.  (4-28)  simplifies  even  further  to  yield  the 
final  result  of  the  discrete  MAP  estimate  of  0  as 
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f 


k  k 

Z  a.r.sinG(r)  +  Z  B.r.cosQ(r) 


j  =  l 


3  3 


irN 


qkb 


j  =  l 


3  3 


F  (Q )  =  0  ,  k  >  3  . 


(4-36) 


For  the  particular  case  when  the  phase  6  is  a  Gaussian  ran- 

2 

dom  variable  with  variance  less  than  0.8  rad  ,  the  density 
function  is  practically  given  by 


f 0  (G )  =  [2TTa^r1/2exp[-6V2o£]  , 


. 2  ,-2, 

0- 


(4-37) 


in  the  interval  (— tt  , tt )  .  Then, 


F  (G )  =  ^[lnf0(0)]  =  -  ~  . 

°0 


By  letting  0  =  0  (r) ,  Eq.  (4-36)  becomes 


6  (r )  = 


qkba  2 
TTN_ 


(4-38) 


k  k 

{-  Z  B.r.cosG(r)  -  Z  a.r.sinG(r)}  .  (4-39) 


j=l 


3  3 


3=1 


3  3 


Equation  (4-39)  is  the  MAP  estimate  of  a  Gaussian  phase  G  in 
terms  of  discrete  observations.  The  MAP  estimate  is, 

No 

therefore,  a  function  of  the  detector  noise  -j-  given  by  Eq. 

2 

(4-3) ,  the  variance  of  the  random  phase  and  the  amplitude 
b  of  the  signal.  Equation  (4-39)  is  of  the  form 


x  =  Ucosx  +  Vsinx 


(4-40) 


and  cannot  be  reduced  any  further.  It  can  be  implemented 
in  the  form  of  a  photon  processor  as  shown  in  Figure  6. 
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GAIN 


Fig  6.  Photon  Count  Processor  for  MAP  Estimation  of  0 

Equation  (4-40)  can,  however,  be  solved  numerically  on  the 
computer.  An  effective  way  is  to  expand  the  sine  and 
cosine  functions  into  series  and  solve  the  resulting  poly¬ 
nomial  iteratively  as  a  predictor-corrector  algorithm.  A 
fifth  order  expansion  (three  terms  in  each  function)  results 
in  reasonably  good  solutions  up  to  0.7  radians,  but  it  gives 
gross  errors  for  higher  phase  angles.  On  the  other  hand,  a 
seventh  order  expansion  (four  terms  in  each  function)  proved 
to  give  very  accurate  results  regardless  of  the  phase  value. 
Thus,  using  four  terms  in  the  expansion  of  the  trigonometric 
functions  of  Eq.  (4-40)  (Ref  12:472),  the  following  poly¬ 
nomial  is  obtained: 
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5040U  +  5040 (V-l ) x  -  2520Ux2  -  840Vx3 
+  210Ux4  +  42Vx5  -  7Ux6  -  Vx7  =  0  (4-41) 

Equation  (4-41)  has  seven  roots.  It  was  solved  using  a 
subroutine  from  the  International  Mathematical  and  Statisti¬ 
cal  Library  (IMSL)  package  available  for  use  with  the  CDC 
6600/Cyber  74  computer  system  at  the  Air  Force  Institute  of 
Technology.  The  real  roots  so  found  that  satisfy  Eq.  (4-40) 
are  the  estimator  solutions.  It  was  found  from  several 
tests  (no  proof  included)  that  the  lowest  root  always  gives 
the  solution  with  the  smallest  error,  and  is  therefore  the 
sought  MAP  estimate.  As  an  example,  the  solution  to 

x  =  11.23108494624  cosx 

-38.93776905131  sinx  ,  (4-42) 

taken  from  an  actual  simulation  run,  is,  using  single  pre¬ 
cision  , 


x  =  0.2740523040524  radians  (4-43) 

with  a  discrepancy  of 

7. 920348821244 E- 9  radians  ,  (4-44) 

which  is  an  extremely  accurate  result. 

For  the  purpose  of  comparison,  it  is  easily  found 
using  Eq.  (4-22)  that  the  MAP  phase  estimate  from  a  continu¬ 
ous  measurement  is  (Refs  11;189  and  15:129) 
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2bu ' 


e  = 


it 
-( 


r(t)cos  ut  dt]cosO 


I 


r(t)sin  tot  dt]sinO^  .  (4-45) 


In  order  to  give  an  idea  of  the  forms  assumed  by  Eq. 
(4-39)  for  a  specified  number  of  observations,  a  few  illus¬ 
trative  examples  are  given  in  Table  I.  Finally,  in  order 
not  to  disrupt  the  continuity  of  the  present  developments, 
verification  of  Eq.  (4-39)  through  computer  simulation  is 
reserved  for  the  last  section  of  this  chapter. 


Maximum  Likelihood  Estimate 

The  ML  estimate  of  a  parameter  6  observed  in  noise  is 
based  on  the  conditional  density  of  the  observations  given  6 
only.  This  is  performed  when  the  statistical  description  of 
theta  is  not  available  (Ref  14:65).  It  is  reasonable  to 
assume  that  0  is  equally  likely  to  occur  in  the  interval 
( — tt  ,  tt  )  not  having  any  other  a  priori  information.  The 
density  of  9  is  therefore  modeled  as 

f0  (0)  =  Yv  '  ^  <  8  £  TT 

=  0  ,  elsewhere  ,  (4-46) 

2 

with  zero  mean  value  and  variance  tt  /3.  From  Eq.  (4-46),  it 

2 

is  found  that  F(0)  =  0  (equivalent  to  a  large  o-  in  Eq. 
(4-37)),  and  the  ML  estimate  is  defined  by  Eq .  (3-33). 
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TABLE  I 


Examples  of  the  Maximum  A  Posteriori  Phase 
Estimator  Form  for  k  Observations 


k 

6(rx  . 

..  rk) 

3 

3  / 2  ~  3 

— r—  (2r, -r0-r-, )  k  cos6(r)  -  4(r0-r-1)k  sin9  (r) 

4  1  2  a  o  —  2  2  a  o  — 

fl 

/2  (r,-r3 ) kQcosO (r ) 

-  /2  (r2-r^ ) kQsin0  (r ) 

3/3  ro, 

4  l2(rrr4)+r2- 

r^-rc+r,]k  cosQ (r) 

J  j  b  O  — 

6 

-  |(r2+r3-r5 

-r, ) k  sin0 (r) 
bo  — 

P 

T 

r 

oo 

[ 

/  r (t) coswtdt] k  cos§ 

-  [  1  r (t) sinotdt] k  sin9 

t  p 

Jo  P 

A 

2 

qkba0 

k 

o 

tiN 

o 

A 

2bo0 

p 

No 
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Equation  (4-36)  then  becomes 


k  k 

Y  a.r.sinO(r)  +  Y  £.r.cos6(r)  =  0 


j  =  l 

for  k  >  3  . 


(4-47] 


(4-47) 

reduces 

to 

r  ^ 

—■ 

!  -  i 

8  ■  r  . 

-1 

i=i 

3  3 

— 

tan 

j 

k 

/ 

ML 

I 

a  .  r  . 

L  j=l 

3  3  J 

,  k  >  3  . 


(4-48) 


Equation  (4-48)  is  the  ML  estimate  of  the  phase  based  on 
discrete  observations.  The  condition  that  k  >  3  is  main¬ 
tained  since  for  k  =  1,  the  tangent  of  0  (r )  is  not  defined 
(zero  divided  by  zero)  and  for  k  =  2  the  summations  in  Eqs. 
(4-29),  (4-30),  and  (4-35)  do  not  equal  zero.  Again  for 
comparison,  as  k  becomes  very  large  and  with  F(0)  =  0,  Eq. 
(4-45)  gives  (Ref  15:129) 


=  tan 


f  r(t)  cos  wt  dt 
-1  _0 _ 

rT 

J  r  (t)  sin  wt  dt 

L  0 


(4-49) 


Equation  (4-49)  is  the  ML  phase  estimate  made  from  a  con¬ 
tinuous  measurement. 

Equation  (4-48)  takes  on  different  forms  for  each  value 
of  k.  In  order  to  visualize  how  it  changes,  a  few  illustra¬ 
tive  examples  are  given  in  Table  II.  These  results  agree 
with  the  findings  of  Wyant  (Ref  16:264),  which  inspired  the 
work  of  this  thesis,  although  he  assumed  shot  noise  instead. 
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TABLE  II 


Examples  of  the  Maximum  Likelihood  Phase 
Estimator  Form  for  k  Observations 


I 

tan6  (r.  ...  r,  ) 

1  k 

3 

✓3  2rrr2_r3 

3  r2”r3 

4 

rl"r3 

r2  r4 

6 

/3  2(rrr4)+r2-r3'r5+r6 

3  r2+r3_rS~r6 

0.3129  (r1-r11)  +  0. 2980  (r2-r10*-r12+r2()) 

0. 0967  (r2+r1Q-r12-r20)  +  0 . 1839 (r3+rg-r13~rig 

20 

+0.2531  (r3-rg-r13  +  r19)  +  0 . 1839  (r4-r8-rig+r18 ) 
•••+0.2531  (rg+r8-r14-r18)  +  0 . 2976  (r^-r^-r^  )  *•* 

...+0-0967(r5-r7-r15+r17) 

+  0.3129  (r  -r, c) 

b  lb 

rT 

J  r  (t)cos  wt  dt 

oo 

C 

,T 

J  r  (t)  sin  tot  dt 

O 
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He  originally  noted  that  k  >  3,  and  proposed  four  inte¬ 
grations.  The  result  corresponds  to  the  entry  for  k  =  4  in 
Table  II. 

Before  proceeding  to  the  algorithm  verification,  it  is 
necessary  to  determine  (or  predict)  how  good  an  estimate 
can  be  obtained  from  Eqs.  (4-36)  and  (4-48).  Therefore, 
the  question  of  performance  remains  to  be  addressed. 

The  conditions  for  which  photon  counting  yields  accu¬ 
rate  estimates  are  to  be  investigated  in  the  next  sections. 

Performance  of  the  Estimators 

One  measure  of  performance  of  an  estimator  is  generally 
given  by  the  variance  of  the  estimate  with  respect  to  the 
true  value.  This  measure  of  performance  is  known  as  the 
Mean-Square  Error  (MSE) ,  and  in  the  case  of  unbiased  esti¬ 
mates  (zero-mean  errors) ,  it  is  the  variance  of  the  error 
itself . 

For  nonlinear  estimates,  the  MS  error  is  not  simple  to 
compute,  and  a  lower  bound  on  the  error  is  sought  instead. 
The  most  widely  used  bound  is  the  one  given  by  the  Craraer- 
Rao  inequality  (Ref  14:66-73).  An  approximate  MS  error 
and  the  Cramer-Rao  (lower)  Bound  (CRB)  will  be  computed  in 
this  paper  and  a  comparison  will  be  made  to  establish  the 
validity  of  the  CRB  bound  to  measure  the  performance  of  the 
joint  estimation  problem  addressed  in  Chapter  VI. 

Mean-Square  Error.  The  MSE  of  an  estimate  is  defined 
as  the  expected  value  of  the  square  error  between  the 
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estimate  0 (r)  and  the  actual  value  of  0  (Refs  14:64,  17:412). 


Let  the  error  be  defined  as 

er  ±  (0-0 (r) )  . 


(4-50) 


The  exact  MSE  is  then  given  by  (Ref  14:56) 


.00  oo 


E[c‘]  = 


j  d0  j  dr(0  -  0 (r) )2fQ^r  (0,r)  .  (4-51) 


—  00  — OO 


The  joint  density  is  given  by  the  product  f r | Q (r | 0 ) f ^ (0 ) . 

For  the  MAP  estimate,  0 (r)  and  the  two  densities  are  given 
by  Eqs .  (4-39),  (4-14),  and  (4-37).  For  the  ML  estimate, 

0 (r )  and  the  densities  are  given  by  Eqs.  (4-48),  (4-14),  and 

(4-46) .  The  direct  computation  of  Eq.  (4-51)  is  difficult 
to  perform.  However,  an  indirect  computation  suggested  by 
Sage  and  Melsa  (Ref  11:189)  is  performed  here  for  discrete 
observations  and  yields  a  result  that  depends  on  the  energy 
of  the  system  for  large  signal-to-noise  ratios.  These  deri¬ 
vations  are  developed  in  the  following  sections. 

The  MAP-MS  Estimation  Error.  The  MAP  estimate  of  6 
when  the  a  priori  density  is  Gaussian  is,  from  Eq.  (4-36), 


~  k  k 

0 (r)  =  -  E  g.r.k  cos§(r)  -  I  a.r.k  sinG (r)  , 


j  =  l 

where,  from  Table  I, 


3  J  o 


.  A  <5kba0 
Ko  ttN 


3  =  1 


3  3  o 


(4-52) 


(4-53) 
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If  the  expression  for  r^  given  by  Eq.  (3-16)  is  substituted 
into  Eq.  (4-53),  then  the  MAP  estimate  becomes,  in  terms  of 
signal  and  noise. 


0  (r)  =  - 


I  6 • s . (0 ) k  cosO (r ) 
j=l  3  3  ° 


Z  6 .n .k  cos0 (r) 

j_l  3  3  0 


k  ~  k 

Z  a jS j  (0  )  kosin0  (r)  -  Z  oun^k^inO  (r)  , 

(4-54) 

where  s^(0)  is  given  by  Eq.  (3-25).  After  some  processing 
of  Eq.  (4-54),  the  estimate  becomes 


a  ■  I  ■  a  a 

— r-  k  Z  {£3.cos0(r)  +  a.sin0(r)} 

qk  o  j=1  3  -  3 


hT  r  /v  2  ~ 

—  k  Z  (B.a.cosO  cos0  (r)  -  g.sin0  cos0  (r) 
ziiq  O  j  =  ^  J  J  D 


+  afcosG  sin0  (r)  -  a^G^sinG  sinG  (r) } 


k  Z  (6.n.cos9(r)  +  a.n.sin0(r)} 
°  j=l  3  3  33 


(4-55) 


The  equalities  of  Eqs.  (4-30)  and  (4-35)  can  now  be  substi- 

k  k 

tuted  into  Eq.  (4-55).  By  setting  Z  a. 6.  =  Z  a.  = 
k  j=l33j=l3^ 

Z  3 ■  =0  and  expanding  the  products  of  sin0 ■ cos0(r )  and 

3  =  3- 

cos0-sin0(r)  into  sums,  Eq.  (4-55)  becomes 


k  Z  {  62  [sin  (G  +  0  (r)  )  +  sin(G  -  6  (r)  )  ] 


4nq  o  ^ j 


-  a ^ [sin  (0  +  0(r))  -  sin (6  -  0  (r) ) ] 


-  k  Z  {g.n.cosG(r)  +  a.n.sinG(r)}  .  (4-56) 

°  j=l  3  3  33 

The  error  e  defined  in  Eq.  (4-50)  is  expected  to  be  small; 
otherwise,  the  estimate  is  meaningless.  To  ins*  e  this 
condition,  the  approximation 


sin (0  -  0 (r) ) 


(4-57) 


can  be  used  in  Eq.  (4-56) .  By  observing  from  Eq.  (4-2?) 
k  2  k  _ 

that  Z  a  7  =  Z  6^  =  0,  Eq.  (4-56)  can  be  transformed 

j=l  3  j=1  J 

to  yield  an  expression  for  the  error;  so,  it  becomes 


"2no  2rrqk 

2  Yk9 bT  Yk 
Ea°0 


Z  {g.n.cos0(r)  -  a.n.sinG(r)}  , 
j  =  l  3  3  3  3 


(4-58) 


where  new  parameters  have  been  defined  as 


A  .v  V  2.-1 
y,  -  [k  Z  a.]  , 

*  j=l.  3 


(4-59) 


which  is  just  a  number  that  depends  on  the  particular  value 


of  k,  and 
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Ea  ~  ,  t1  <  t  <  tx+T  ,  (4-60) 

which  is  the  ac  component  of  the  energy  in  the  interval 
( ti , ti+T)  of  the  signal  of  Eq.  (2-9). 

A 

The  first  check  to  be  made  on  the  estimate  0 (r)  using 
Eq.  (4-58)  is  that  for  bias  (Ref  17  :  404)  .  Thus,  by  recall¬ 
ing  from  Eq.  (4-7)  that  E [ n j ]  =  0,  the  conditional  expecta¬ 
tion  of  the  error  given  0  is 

7T2N 

E[£r|0]  =  - fYkE[0U)|e]  .  (4-61) 

Eaa0 

A 

Solving  for  E[0(r)|0]  and  noticing  that  E [ 0 | 0 ]  =  0,  Eq. 

(4-61)  gives 

E  [ 0  (r )  |  0  ]  =  0[1  +  TT2NoYk/Eaa2]"1  .  (4-62) 

A 

Thus,  the  conditional  expected  value  of  0  (r)  is  a  scaled 
version  of  the  true  value  of  0.  The  estimate  approaches 
the  true  value  as  the  signal  to  noise  energy  ratio  in 
(ti,ti+T),  Eg/No,  is  made  very  large,  which  is  expected  when 
the  measurement  is  made  in  the  absence  of  noise.  Thus,  in 
the  limiting  case,  and  only  in  that  case, 

lira  E[§  (r) 1 0]  =  0  ,  (4-63) 

E  /N 
a  o 

and  the  MAP  phase  estimate  is  asymptotically  unbiased  (Ref 
17:404) . 
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The  second  computation  to  be  made  from  c is  that  of 

the  MS  error  between  0 (r)  and  0.  The  MS  error  is  defined 

2 

as  the  expected  value  of  r-r.  The  first  order  statistics 
needed  to  asymptotically  characterize  the  MSE  are 


E[0] 
E  [0  (r)  ] 
E[er] 


0 

E[E[0  (r)  |  0]  ]  =  0 

E[0]  -  E[0  (r) ]  =  0  . 


(4-64) 

(4-65) 

(4-66) 


In  such  a  case,  large  E  /N  ,  the  MS  error  is  both  the  vari- 

a  o 

ance  of  0 (r)  and  the  variance  of  er  which  will  be  defined 
asymptotically  as 


V  (s  )  £  lim  V  ( e  ) 

a  r  E  /N  ar  r 

a  o 


=  lim  E  [E‘]  . 

E  /N  -*•<» 
a  o 


(4-67) 

From  Eg.  (4-58)  the  second  order  moment  of  er  is  given  by 


E[  A  =  [tt2N  Y. /E  a?]  2  E[(9(r))2] 

L  O  A.  d  C 

+  [2irqkYk/bT]  2 

E[{  Z  (8.n.cos0(r)  +  a . n . sin§ (r ) ) } z ]  ,  (4-68) 

j=l  3  3  33 


where  the  identity  L[n^]  =  0  was  used, 
tion  term  in  Eg.  (4-68)  reduces  to 


The  second  expecta- 


N  T  k  7 

-V  1  aj 

2gzk  j=l  3 


t 
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where  use  of  the  identities  of  Eqs.  (4-8)  ,  (4-29)  and  (4-30) 
was  made.  Finally,  by  using  Eqs.  (4-59)  and  (4-60),  Eq. 
(4-68)  becomes 


EUJ)  =  [TT2NoYk/Eaa^]2E[(0(r))2] 


+  U  Vk/Ea  * 


(4-69) 


Equation  (4-69)  is  the  MS  error  of  the  phase  estimate.  The 

A  2 

difficulty  with  this  expression  is  that  E[(0(r))  ]  must  be 

evaluated  from  the  density  function  of  Eq.  (4-14)  (Ref  14:64) 
However,  if  the  ratio  Ea/NQ  is  made  very  large,  the  error 


variance  becomes,  from  Eq.  (4-67) 


Va(er}  = 


tt2N  y, 
o'k 


(4-70) 


It  will  be  beneficial  to  inquire  what  the  effect  of  y^  in 
Eq.  (4-70)  is  as  k  varies  from  k  =  3  to  the  limit  where  it 
becomes  very  large.  A  plot  of  y^  versus  k  for  large  signal- 
to-noise  ratio  (SNR)  is  shown  in  Figure  7  as  a  normalized 


variance 


VEr,Ea 


2-1 

A  plot  of  y,k;  that  is,  [  Z  a.]  ,  versus  k  is  shown  in 

*  j=l  3 

Figure  8  to  show  how  fast  the  (normalized)  variance,  y^, 
approaches  linearity  as  k  increases.  From  Figure  7  it  is 
observed  that  y^  tends  to  stabilize  at  around  0.051  as  k 
becomes  large.  This  corresponds  to  the  asymptotic  slope  of 
the  curve  in  Figure  8.  Thus,  equation  (4-70)  becomes 
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Fig  8.  Relationship  between  y,  and  I<  as  a  Product 
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1 . i .m. 

k  ->  oo 


V 
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(e  ) 
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0 . 051tt  N 


o 


(4-71) 


The  MS  error  of  a  phase  estimate  for  large  signal-to-noise 
energy  ratio  in  (t^t^+T)  is  just  a  multiple  of  the  inverse 
SNR  scaled  by  a  factor  that  depends  on  the  number  of  the 
observation  interval  subdivisions,  k.  A  few  examples  to 
illustrate  Eq.  (4-70)  are  given  in  Table  III. 

The  ML-MS  Estimation  Error.  The  same  procedure  used 
in  the  previous  section  to  determine  the  MAP-MS  error  can  be 
used  for  the  ML  case.  Using  Eq.  (4-48)  as  a  starting  point 
and  substituting  for  r ^  given  by  Eq.  (3-16),  it  is  found  that 
the  ML  estimate  is  unconditionally  unbiased;  that  is, 


E  [6  (r)  I  0 ]  =6 


(4-72) 


This  is  easy  to  verify  by  observing  that  for  ML,  of  =  °°  in 

o 

the  MAP  equations.  Thus,  Eq.  (4-69)  reduces  to  Eq.  (4-72) 
without  any  SNR  restrictions.  Completion  of  the  procedure 
yields  an  unconditional  error  variance  given  by 


v,r(tr>  ’ 


tt2N  Yv 
o  k 


(4-73) 


Again  this  is  easily  verified  by  letting  =  <»  in  Eq. 
(4-69).  A  different  procedure,  however,  carried  out  by 
Raemer  (Ref  9:263,458)  to  compute  directly  the  Root  Mean- 
Square  (RMS)  error  defined  as 
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TABLE  III 


Examples  of  the  MAP  Phase  Estimator 
Error  Variance  for  Large  SNR 


k 

v  <e  ) 

a '  x' 

3 

.2  N  /2  N  /2 

4  7T  o'  -(  o' 

27  E  X'462  E 

a  a 

4 

2  N  /2  N  /2 

■no  _  ^  2337  ° 

8  E  E 

a  a 

5 

0  2  N  /2  N  /2 

2r  o  i  1427  ° 

17.27  E  1,1  E 

a  a 

6 

2  N  /2  N  /2 

v  -f—  - 

a.  3 

15 

02  N  /2  N  / 2 

2tt  o  o 

19.45  E  “  1.0147 

a  a 

20 

~  2  N  /2  N  /2 

19.57  E  =  1-0082  E 

3  3 

00 

N  /2  N  /2 

o  _  o 

E  E 

a  a 
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[V  (c  )] 
ar  r 


1/2 


k  =  00  , 


(4-74) 


where  0  (uo)  is  as  in  Eq.  (4-49)  did  show  a  SNR  dependence  of 
the  ML  estimate.  Dy  expanding  the  arc  tangent  in  Eq.  (4-49), 
Raemer  obtained  similar  results  with  the  condition  of  a 
large  Ea/NQ  ratio.  Specifically, 


4  p . BnT 
i 


N_o 

2E. 


(4-75) 


where  and  Bn  were  defined  parameters  of  the  signal  and 
noise  treated  as  narrowband  processes.  Therefore,  the  MAP 
and  ML  error  variance  are  the  same  and  approach  zero  when 
the  signal  is  strong  relative  to  the  receiver  noise.  Thus, 
Table  III  is  equally  an  ML  comparison  of  the  error  variance 
in  terms  of  the  observations  k.  To  show  graphically  the 
relative  improvement  of  entries  in  Table  III,  define  in 
logarithmic  terms 


I  (dB)  -  10  log 


Var(er>  k=3 


V  (e  ) 
ar  r 


=  10  log  (Y3/Yk) 


(4-76) 


The  relative  improvement  in  decibels  (dB)  of  the  error  er 
with  respect  to  the  error  in  three  subintervals  is  shown 
in  Figure  9. 

The  important  conclusion  of  this  section  is  that  in 
the  presence  of  thermal  (receiver)  noise,  the  error  in  the 
phase  estimate  is  strictly  a  function  of  the  signal-to-noise 
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ratio.  The  ideal  procedure  is  to  perform  a  current  measure¬ 
ment  (k  =  °°)  ,  for  which  the  error  is  the  inverse  SNR  multi¬ 
plied  by  a  factor  of  1.  The  procedure  is  suboptimal  if 
discrete  photon  counts  are  made  (k  <  «>)  ,  for  which  case  the 
same  error  is  scaled  by  an  appropriate  factor  determined  by 
k.  Referring  to  the  curve  in  Figure  9,  the  continuous  cur¬ 
rent-based  estimate  offers  an  improvement  of  1.65  dB  over 
the  discrete  photon  count  based  on  only  k  =  3  integrations. 
Notice,  however,  that  if  the  hypothetical  CCD  detectors  can 
provide  15  or  more  integrations  over  T,  then  the  improve¬ 
ment  of  a  continuous  measurement  over  photon  counting  is 
only  0.064  dB.  Thus,  provided  the  SNR  restrictions  are  met, 
the  photon  counting  algorithms  offer  a  viable  alternative 
for  phase  estimation  when  the  signal  levels  are  low. 

Considering  the  extremes,  when  the  SNR  is  low,  neither 
a  current  measurement  nor  a  photon  count  will  yield  an  accu¬ 
rate  phase  estimate.  On  the  other  hand,  in  the  absence  of 
noise,  the  error  is  zero  (E^/N  =  <*>)  as  expected  from  a 

perfect  measurement.  This  is  easily  proven  by  substituting 
r j  =  Sj(0),  given  by  Eq.  (4-23),  back  into  Eqs.  (4-36) 

(with  F(0)  =  0)  or  (4-48).  With  the  help  of  Eqs.  (4-29), 
(4-30) ,  and  (4-35)  the  result  obtained  is  6 (r)  =  0  as 
expected.  This  was  also  verified  computationally. 

The  level  of  SNR  to  assure  with  a  given  certainty  that 
an  estimate  is  within  a  specified  phase  range  has  been 
addressed  in  the  literature  and  will  not  be  repeated  here 
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(Ref  9:266).  Instead,  a  general  observation  will  be  made 
from  the  results  of  the  computer  simulation  discussed  at 
the  end  of  this  chapter. 

In  order  to  complete  the  error  analysis,  it  is  desir¬ 
able  to  inquire  what  the  minimum  MS  error  is,  and  the  con¬ 
ditions  under  which  it  is  attained.  A  good  description  is 
given  by  the  Cramer-Rao  bound  on  the  minimum  MSE  although 
other  bounds  are  also  available  (Ref  14:71).  The  Cramer- 
Rao  approach  will  be  considered  in  the  next  section. 

Cramer-Rao  Lower  Bound  on  the  Minimum  MS  Error.  The 
lower  limit  on  the  value  that  the  minimum  MS  error  (MMSE) 
can  have  can  be  calculated  without  actually  having  to  know 
the  estimate.  For  the  case  of  an  unbiased  estimate,  the 
CRB  bound  is  given  by  (Ref  14:72) 

CRB  =  {E[  (^  Inf  Q  (r,0))  2]  J"1 

82  -1 

=  - { E [  — — y  Inf  ( r , 0 ) ] }  .  (4-77) 

ae2  *>e  - 

For  the  case  of  a  biased  non-random  parameter,  a  bound 
expression  is  given  by  Van  Trees  (Ref  14:147). 

In  order  to  use  Eq.  (4-77)  in  the  problem  at  hand,  it 
is  necessary  to  work  with  the  large  SNR  restriction,  for 

-'N 

which  the  estimate  9  (r)  was  shown  to  be  asymptotically 
unbiased.  Thus,  the  bound  s  given  by 

V  (c)  £  CRB  ,  E  /N  ►  ~  .  (4-78) 

a  IT  ci  O 
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Before  letting  the  SNR  become  very  large  let 


var(cr)  >  CRB  =  ~{E[|o'  lnfr|G(-i0)  +  Inf^O)]}-1 

(4-79) 

For  the  observation  in  white  Gaussian  noise,  it  is  conven¬ 
ient  to  specialize  Eq.  (4-79)  further  before  considering 
the  phase  problem.  From  Eq.  (4-17) 


36  lnfr|0 


( r  |  9 )  = 


2-3.  Ji  £  r  JL  s  (q  \ 

NqT  rj  30  Sj(0) 


-  T^r  7,  V9’  A  V9)  ' 

o  3=1  J  J 

(4-80) 

The  second  derivative  of  Eq.  (4-80)  with  respect  to  0  is 


30 


2  lnfr|0(r-l0)  = 


0  2,  k  l 
ZaJs  z  )  rr 
NT  .  .  ) lrj 

o  3  =  1  l  J 


-  3^  (0)  3 


_3J 

36 


i 


-m' 


(4-81) 


The  CRB  bound  then  becomes 


CRB  = 


£  E 

I  NoT  j  =  l 

.  (?V0)V 

V  36  / 


(r.  -  Sj (0) ) 


+  E 


30 


2  Sj(6) 


"I  i  -1 


30 


2  lnfQ 


(4-82) 


By  observing  that 
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and 
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(4-85) 


By  use  of  Eq .  (4-85)  in  Eq.  (4-83),  the  bound  becomes 


hq2k  bT  2 


CRB  *  ji^f  <j^>‘  E  i“jsin0  +  V030’2  -  ei|t  ftoijj  1 


(~)2  Y.  [a2sin20  +  62cos20  +  a.8.sin20) 
|  No  j=l  3  3  33 


-  E[g^  F(0)]j  ,  k  >  3  . 


(4-86) 


Invoking  the  identities  of  Eqs.  (4-29)  and  (4-30)  one  more 
time,  Eq.  (4-86)  can  be  simplified  to 


Va<er] 


kb  T 


I-2 


4 it  Nq/2  j  =  l 


k  2  3  l“] 

£  “2  -  F(e)1( 


E,  a  | -1 

IT—  -  E[re  F(0,,(  '  k  *  3  ' 


•"‘Vk 


for  E  /N  ■+  00  .  (4-87) 

a  o 

Equation  (4-87)  is  the  lower  bound  on  the  MS  error  in  the 
estimate  of  0  given  by  the  Cramer-Rao  inequality.  It  is  a 
function  of  the  ac  component  of  the  signal  energy  E  .  the 
noise  energy  Nq  and  the  observation  subinterval  k  repre¬ 
sented  by  .  It  is  also  a  function  of  the  statistical 
description  of  6  represented  by  the  density  fg(0).  In  the 
particular  case  where  the  phase  is  modeled  as  a  Gaussian 
random  variable  with  the  density  given  by  Eq.  (4-37),  then 
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(4-88) 


EIw  F<°>>  '  -  -T  ' 

°6 

The  bound  on  the  error  variance  becomes 


V  (e  ) 
a  '  r 1 


°e 


-l 


2kl  2 

71  Vk°e 


'2"oVe  +  Ea°e  ' 


k  >  3 


for  E  /N  °°  . 
a  o 


(4-89) 


In  the  limiting  case  as  k  becomes  very  large,  1/y,  ^  19.6 

(see  p.  44) ,  so 


l.i.m.  V  (e  )  > 

,  a.  I. 

k  -*•  <» 


2E  . 
a  1 
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°  aeJ 


-1 


2 
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N  +  2E  o‘ 
o  a  0 


for  Ea/No 


(4-90) 


Equation  (4-90)  gives  the  MAP  lower  bound  on  the  error  var¬ 
iance  based  on  a  continuous  measurement.  Finally,  when  the 
large  SNR  condition  is  used,  Eq.  (4-90)  becomes 


N 

lim  {l.i.m.  V  (e  )  }  i 

E  /N  ->»  k  ®  a  r 
n  o 


(4-91) 


This  is  precisely  the  result  given  by  Eq.  (4-71) .  There¬ 
fore,  as  the  term  E  /N  grows  unbounded,  the  true  error 

ci  O 
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variance  approaches  the  CRB  with  equality;  the  MAP  estimates 
0 (r )  and  6  are  asymptotically  efficient  (Ref  14:276).  Thus, 
the  Cramer-Rao  bounds  are  also  given  in  Table  III.  The 
same  observation  is  made  for  the  ML-CRB  on  the  error  vari¬ 
ance.  For  the  ML-CRB,  from  Eq.  (4-46), 

E[^g-  F(0)  )  =  0  ,  (4-92) 

and  Eq.  (4-77)  reduces  to  (Ref  14:66) 

a2  -1 

CRB  =  —  { E  [  — — Inf  lfl(rjO)]}  .  (4-93) 

36 

Therefore  the  bound  for  the  ML  error  is,  from  Eq.  (4-89)  as 
2 

oQ  gets  large, 

*2Vk 

Va(er}  ~  . E  '  k  -  3  '  (4_94) 

3. 

or  from  Eq.  (4-90)  , 

N 

l.i.m.  Va(er)  >  j~  .  (4-95) 

k  •*  “  ^a 

From  Eq.  (4-90)  it  is  observed  that  when  the  signal  is 
,  2 

weak,  the  error  is  limited  by  the  variance  o0.  But  this  is 
not  the  same  result  expressed  by  Eq.  (4-69) ,  and  therefore 
the  Cramer-Rao  bound  becomes  meaningless  in  such  a  case. 

In  conclusion,  the  analysis  made  in  the  Gaussian  noise 
context  indicates  that  under  thermal  noise  limited  condi¬ 
tions,  measurement  of  a  phasefront  cannot  be  made  accurately 
except  for  large  signal-to-noise  ratios.  For  low  light 
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levels  conforming  to  this  restriction,  the  photon  counting 
technique  is  quite  acceptable,  but  when  the  signal  level 
permits,  a  current  measurement  is  most  appropriate.  At 
this  point,  a  verification  of  the  algorithms  derived  in 
Eqs.  (4-39)  and  (4-48),  as  well  as  a  performance  evaluation, 
is  overdue.  This  was  done  in  a  computer  simulation,  and  is 
presented  in  the  following  section. 


Verification  of  the  MAP  and  ML  Estimators 

A  simulation  program  to  test  Eqs.  (4-39)  and  (4-48)  was 
written  in  FORTRAN  and  run  on  the  CDC  6600  computer  system. 
After  the  parameters  for  the  simulation  have  been  input, 
the  program  generates  the  signal  counts  s^(6),  computed 
from  the  true  phase  input  and  Eq.  (3-24),  and  adds  white 
Gaussian  noise  counts  n ^ ,  computed  from 

P.  NT  1/2 

n.  =  -Z  [-£-]  ,  (4-96) 

j  q  2k 


where  P .  is  a  number  from  a  zero-mean,  unit-variance, 

3 

Gaussian  random  number  sequence,  generated  using  a  subrou¬ 
tine  from  the  IMSL  library.  The  program  then  estimates  the 
phase  based  on  the  noise -corrupted  measurements  r_.  using 
Eqs.  (4-39),  (4-41)  and  (4-48)  as  needed. 

The  parameters  chosen  for  this  simulation  are  a  25 
kilohertz  modulation  frequency,  a  0.04  millisecond  observa¬ 
tion  period  (required  by  the  choice  of  N  =  1  in  Eq.  (3-22) ) 
and  a  signal  level  in  the  order  of  1  microamp.  Although 
the  algorithm  is  independent  of  signal  level  (also  verified 
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by  using  10  milliamps)  this  number  was  entered  to  simulate 
the  output  of  a  detector  under  low  signal  levels.  The  num-- 

O 

ber  of  counts  r^  thus  generated  are  of  the  order  of  10 

photoelectrons  and  below.  In  the  0.04  millisecond  observa- 

13 

tion  time  used,  this  corresponds  to  less  than  10  photo¬ 
electrons  per  second,  a  rate  which  is  in  the  category  of 
low  photon  coherence.  The  simulation  results  are  contained 
in  Figures  10-19. 

Simulation  Results  for  the  MAP  Estimator.  The  MAP 

estimator  simulation  is  a  series  of  tests  performed  on  Eq. 

(4-39)  to  verify  its  performance  as  an  estimator.  In  these 

tests,  Eq.  (4-39)  was  used  to  estimate  a  known  phase  given 

different  signal-to-noise  ratios  following  the  procedure 

outlined  in  the  previous  section.  The  tests  were  then 

repeated  to  estimate  different  phase  angles  at  a  fixed  SNR. 

Several  such  tests  were  performed  and  the  results  presented 

here  are  typical.  The  first  set  is  shown  in  Figure  10. 

These  curves  were  obtained  using  a  phase  variance  of  0.8 
2 

rad  ,  chosen  as  the  safe  maximum  deviation  for  which  Eq. 
(4-38)  can  be  used  to  obtain  Eq.  (4-39);  and  using  a  test 
phase  of  0.5234  radians  (30°),  which  seemed  a  logical  choice. 
The  five  curves  shown  are  the  phase  estimates  plotted  as  a 
function  of  the  algorithmic  form  used  (determined  by  k)  for 
the  indicated  signal-to-noise  ratios  (10  dB  to  30  dB) .  The 
same  noise  counts  were  used  in  each  curve. 

From  this  and  other  tests  performed  (by  changing  the 
noise  seed  to  generate  P^),  it  is  observed  that  the  estimate 
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WHITE  OfiUSSIRN  NOISE 


Fig  10.  MAP  Simulation  Results  -  Single  Phase  at  Different  SNR' 


deviations  are  larger  for  lower  values  of  k,  but  improve  as 
the  SNR  increases.  Thus,  for  large  SNR's,  consider  that  to 
be  50  (17  dB)  and  above,  the  benefit  of  the  simpler  algo¬ 
rithm  forms  in  Table  I  can  be  used  to  an  advantage  (less 
computation  time) .  The  same  is  true  for  the  ML  approach 
(see  Table  II) .  The  variance  of  these  estimates  is  shown 
in  Figure  11. 

Figure  12  is  a  plot  of  the  estimates  of  five  phase 
angles  obtained  with  a  SNR  of  17  dB.  They  range  from  0.2 
radians  (11.5°)  to  1.4  radians  (80°).  No  significant  dif¬ 
ferences  in  the  pattern  are  observed.  A  small  discrepancy, 
however,  had  to  be  fixed  in  this  particular  plot;  the  phase 
estimate  0(r4)  for  0  =  1.4  radians  made  by  using  four  sub¬ 
intervals  had  a  negative  sign.  This  occurred  both  in  the 
MAP  and  the  ML  test.  However,  it  cannot  be  inferred  that 
the  estimator  is  more  sensitive  as  the  phase  becomes  larger. 
The  variance  of  these  estimates  is  shown  in  Figure  13.  Both 
Figures  11  and  13  seem  to  nfirm  the  theoretical  prediction 
of  Figure  7  :  improve^  performance  as  k  (and  SNR)  -> 

A  test  to  check  the  effect  of  a  smaller  a  priori  phase 
2 

variance  was  also  run  and  is  shown  in  Figure  14.  The 

2 

variance  was  chosen  to  be  0.274  rad  (the  square  of  the 
test  phase) .  A  slight  improvement  was  obtained  with  respect 
to  the  curves  of  Figure  10.  The  comparison  can  be  made  more 
easily  by  looking  at  the  estimates  variance  shown  in  Figure 
15.  Here,  the  improvement  is  more  obvious  for  the  lower 
values  of  k,  and  in  particular  for  the  10  dB  curve;  that  is, 
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for  the  estimates  obtained  with  the  smaller  SNR  used.  This 
seems  to  confirm  the  discussion  pertaining  to  Eq .  (4-90), 

that  estimates  from  noisy  measurements  are  more  dependent 
on  the  a  priori  information,  a  useful  characteristic  of  the  MAP 
estimator.  Similar  tests  were  also  run  to  verify  the  ML 
algorithms,  and  are  discussed  in  the  following  section. 

Simulation  Results  for  the  ML  Estimator.  The  ML  esti¬ 
mator  simulation  is  a  series  of  tests  performed  on  Eq. 

(4-48),  similar  to  the  ones  performed  on  the  MAP  estimator. 

The  main  difference  lies  in  the  straightforward  computation 
given  by  the  inverse  tangent  form  of  Eq.  (4-48).  Here,  no 
a  priori  information  (phase  variance)  is  used.  The  results, 
shown  in  Figures  16-19,  are  surprisingly  close  to  the  MAP 
results  already  discussed,  especially  as  k  increases.  Thus, 
the  ML  approach  can  be  used  to  advantage  when  the  signal- 
to-noise  ratios  are  high  enough  to  perform  nearly  as  well  as 
the  MAP  estimator.  From  these  results,  Figures  10-19,  it 
can  be  seen  that  high  enough  may  be  10  dB  and  higher.  Thus, 
given  SNR's  better  than  10  dB,  the  ML  algorithm  of  Eq.  (4-48) 
seems  to  be  a  good  estimator.  Its  use,  when  warranted,  has 
the  advantage  of  avoiding  solving  equations  of  the  form  of 
Eq.  (4-41)  which  are  time  consuming.  The  MAP  estimator, 
on  the  other  hand,  is  more  useful  under  noisy  conditions, 
where  the  estimator  weights  more  heavily  on  the  a  priori 
information . 

In  the  analysis  presented  in  this  chapter,  the  ultimate 
noise  limited  condition  was  assumed:  that  of  the  detector. 
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Fig  16.  ML  Simulation  Results  -  Single  Phase  at  Different  SNR' 


Simulation  Variance  -  Single  Phase  at  Different  SNR' 


6NR=60 


» 
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Fig  18.  ML  Simulation  Results  -  Multiple  Phase  Angles  at  a  Fixed  SNR 


71 


Fig  19.  ML  Simulation  Variance  -  Multiple  Fhase  Angles  at  a  Fixed  SNR 


In  reality,  and  particularly  for  low  light  levels,  when  the 
signal- to-thermal  noise  energy  ratio  is  indeed  very  large 
as  may  be  the  case  of  the  CC  devices,  then  the  effects  of 
the  signal  shot  noise  predominate.  Therefore,  it  is  neces¬ 
sary  to  consider  the  predominant  shot  noise  case  in  order 
to  have  a  more  complete  analysis.  When  the  phase  0  is  the 
only  random  parameter  in  Eq.  (2-9)  and  consequently  in  Eq. 
(4-23),  the  problem  can  be  treated  in  the  Poisson  context. 
An  explicit  use  of  the  Poisson  statistical  description  of 
shot  noise,  ignored  by  Wyant,  will  be  the  basis  of  the 
developments  of  the  next  chapter. 


V  Phase  Estimation  in  Shot  Noise 


When  detector  thermal  noise  is  negligible,  the  process 
is  signal  shot  noise  limited.  For  a  multimode  field  (Ref  2: 
87-94,  212),  the  process  is  well  modeled  as  having  condi¬ 
tionally  Poisson  statistics.  When  the  field  is  single  mode, 
then  it  is  governed  by  Laquerre  (Ref  2:304)  statistics  and 
will  not  be  considered  here.  The  analysis  is  identical  to 
the  one  made  in  Chapter  IV  after  Eqs.  (3-32)  and  (3-33). 

For  multimode  detection,  the  counts  r.  in  the  observation 

3 

interval  (t^,t^+T)  are  independent  and  Poisson  distributed 
when  conditioned  on  the  field  intensity  (Ref  2:295).  The 
density  function  of  the  observables  in  (t^,tj+i)  conditioned 
on  6  is  the  probability  density  of  the  events  (photons 
received) ,  assuming  that  all  the  events  equal  the  observa¬ 
tions  (photoelectrons  produced)  in  each  observation  sub¬ 
interval  x.  Therefore,  the  density  function  of  the  shot 
process  is 


P[Nx 


(5-1) 


where 


t .  +T 

,  -  /’ 


E[r  (t)  ] 


dt 


(5-2) 


The  observation  r(t)  has  the  form 


r(t)  =  s  ( t ,  0 )  +  q(Agn(t)  +  +  ^d(t))  ,  (5-3) 
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where  the  noise  count  (photon)  rates  have  the  following 

statistics:  (1)  the  signal  shot  noise  A  (t)  is  a  non- 

sn 

stationary,  zero-mean  process  with  covariance  A(t)6(t-t'); 

(2)  the  background  is  assumed  constant  and  infinite  with 

correlation  function  N  ,  6 (r  -r  ')6(t-t')  and  a  stationary 
count  rate  A^  =  N0bB0D0n/hf  ,  where  is  the  spectral 

background  noise  strength,  Bo  is  the  optical  bandwidth  and 
Dg  is  the  number  of  spatial  modes  (Ref  2:212,213,298);  and 

(3)  the  stationary  detector  dark  current  has  dc  content  Ad 
and  covariance  A^i(t-t'). 


Maximum  A  Posteriori  and  Maximum  Likelihood  Estimates 

With  the  above  conditions  established,  the  MAP  and  ML 
phase  estimates  can  be  found  using  Eq.  (5-1)  as  the  starting 
point.  Thus,  Eq.  (5-1)  becomes 


(s. (0)  +  (Ab+Ad)T/k)  3 


exp [- (s j (9 )  +  (Ab+Ad)T/k)]  .  (5-4) 


In  a  manner  analogous  to  Eq.  (4-14),  the  observation  vector 
r  has  a  conditional  density  function  given  by 


k 
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Unlike  Eq.  (4-14),  Eq.  (5-5)  contains  a  constant  noise  term 
(Xfo+A^)T/k  that  represents  deviation  from  ideal  behavior. 
Use  of  Eq.  (5-5)  in  Eq.  (3-32)  yields  the  MAP  estimator 
equation 

k  r.  3s.  (G)  /'T 

4,  S . ( 0 )  +  U.+A  )T/k  30  q  /  ¥0  s(t'6)dt 

J  i-  J  DU  Vq 

+  F (0 )  =  0  ,  (5-6) 


where  the  substitution 


k  3s.  (9) 

E  - i- 

3=1 


30 


i 

j=i  4 


3  s (t, 0) 

30  q 


dt  = 


■  kf 

O 


3s  ( t , 0 ) 


30 


dt 


(5-7) 


and  the  definition  of  Eq.  (4-26)  were  used. 

Use  of  Eqs .  (3-25) - (3-27)  and  (4-24)  in  Eq.  (5-6)  with 

rT  3 

/  tttt  s(t,Q)dt  =  0,  where  s(t,0)  is  given  by  Eq.  (2-9), 

O  30 

yields  the  MAP  estimate 


k  [a . sinO (r ) +6 . cosO (r ) ] k/2r 

L  r  .  - — 4 - 7^ - 

j  =  l  -1  [a  .  cos§  (r) -3  .  sinO  (r )  ]  k/2r+ (a+i,  +i  . ) /b 

J  ~  J  ~  b  d 

-  F ( 9 )  =  0  ,  (5-8) 

where  the  substitution  (Ref  2:113)  (i^+i^)  =  was 

made  to  convert  constant  photon  count  rates  to  ac  noise 
current.  The  ratio 


Y 


A  b 
a  +  ib 


(5-9) 
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is  the  fringe  visibility  and  depends  on  the  radiance  of  the 


source.  The  dark  current  i^  is  usually  negligible. 

Equation  (5-8)  is  the  MAP  discrete  estimator  of  phase 
measured  under  Poisson  shot  noise  conditions.  It  can  be 
reduced  further  under  special  conditions,  but  before 
specializing  its  results  and  following  the  practice  of 
previous  sections,  Eqs.  (5-6)  and  (5-8)  can  be  put  in  the 
form  of  continuous  waveform  equations.  Thus,  in  the  limit 
as  k  becomes  very  large,  Eq.  (5-6)  becomes  (Ref  2:298) 

i  fT  _ -Lit) _  3s(t,0)  _  1  fT  3s  (t,0J 

q  /  s(t,0)  +  (i,+i, )  36  q  /  3e 

J0  o  a  jQ 

+  F(G)  =  0  ,  (5-10) 


and  the  MAP  phase  estimator  of  Eq.  (5-8)  becomes 


n 

if 


r  (t) cos (wt+9 ) 


lQ  sin(wt+0)  +  (a+i^+i^) /b 


dt  -  F (0)  =  0  .  (5-11) 


The  first  case  to  be  considered  is  when  the  fringe 
visibility  is  approximately  unity.  This  is  more  or  less 
what  is  found  from  the  output  of  the  shearing  interferometer, 
as  can  be  determined  from  Eqs.  (2-10)  and  (2-11).  Thus, 
from  these  two  equations,  a/b  -  0.711;  so,  (a+i^)  ~  b,  con¬ 
sidering  the  presence  of  a  small  background  current.  Thus, 
by  setting  (a  +  i^fi^)/b  =  .1,  Eq.  (5-8)  can  be  transformed 
using  simple  trigonometric  identities  (Ref  11:225-241)  into 
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E  r.{tan[sin  ^  ( a  .  /Va"?  +  (5?)  -  y  (r)  ] 
j  =  1  D  3  3  3  ~ 

-1 

+  2n/ [ (a j sinO  (r )  +  3  cosO (r ) ) /k] } 

-  F  ( 0  )  =  0 

for  y  =  id  -  1 

and  Y  >>  id  »  (5-12) 

where  the  substitution 

tan  =  sin'1  (cijAaj  +  3?)  (5-13) 

was  made  to  avoid  dividing  by  zero  as  some  B^'s  are  found  to 
be.  Equation  (5-12)  is  therefore  the  discrete  MAP  phase 
estimator  given  Poisson  shot  noise  conditions  and  equal 
amplitude  dc  and  ac  components.  It  is  noteworthy  to  observe 
that  if  k  is  large,  the  second  term  in  Eq .  (5-12)  can  be 
neglected  so  the  estimate  becomes  approximately 

k  - 

E  r. cot  [sin  ^(ot./Vct?  +  8!p  -  9  (r)  ] 

j = i  ^  a 

-  F  ( 0 )  =  0  ,  k  -*  large  .  (5-14) 


Finally  as  k  becomes  very  large,  from  Eq.  (5-11),  the  con¬ 
tinuous  measurement  estimate  is  found  to  be 
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0 


T 

1 


r  (t) cot ( 


:o  t  +  0 


dt  -  q  F  (0 ) 


for 


Y  +  i 


1  /  Y  >>  i 


d  * 


(5-15) 


The  second  case  to  be  considered  is  when  (a+i^+i^)  =  0 ,  dc 
components  blocked  (capacitively  perhaps).  For  this  case, 
Eq.  (5-8)  becomes 


k  a.sinO(r)  +  5.  cost?  (r) 

Z  r  .  -3 - -- - 2 - 

,  *1  '*•  />. 
j  =  l  a.cos-J  (r)  -  (S.sinC  (r) 


-  F(0)  =  0  . 


(5-16! 


Equation  (5-1G)  can  be  also  transformed  using  simple  trigo¬ 
nometric  identities  into 


Z  r.  tan  [sin  1  (6  •  +  B^)  +  0  (r)  ] 

j  =  l  3  D  J  J 

-  F ( 0 )  =  0  ,  dc  blocked.  (5-17) 


This  is  the  discrete  HYP  phase  estimator  given  Poisson  shot 
noise  conditions  and  dc  components  (signal  and  noise) 
blocked.  In  the  limit  as  k  becomes  large,  the  continuous 
measurement  estimators  found  to  be  from  Eq.  (5-11) 

T 

r(t)cot(o;t  +  6)dt  -  q  F(0)  =  0  , 

dc  blocked  .  (5-18) 


The  last  case  to  be  considered  is  for  a  very  low  fringe 
visibility  such  that  (a+i^)  >>  b.  This  condition  is  given 
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under  a  strong  infinite  background  competing  with  the  source 
(target),  and  is  referred  to  as  background  limited.  Thus, 
for  background  limited  conditions, 

(a  +  i^  +  i^)  /b  >>  [ajCOsG(r)  -  6 ^  sinG  (r )  ]  k/2ii 
and  Eq.  (5-8)  becomes 


£  r .  [a . sinC  (r )  +  p.cosO(r)] 
j=l  3  3  3 


2 if  (a  +  ib  +  id) 
kb 


F  ( Q )  =  0  . 


(5-19) 


Equation  (5-19)  is  of  the  same  form  as  Eq.  (4-36)  with 

the  term  2tt  (a+ib+.>a)/kb  in  place  of  i7No/gkb,  Thus,  an  equiv- 

alent  noise  can  be  defined  as  N  —  2q(a+i.  +i,).  The  phase 

o  d  d 

estimator  in  multimode  shot  noise  under  background  limited 
conditions  is  thus  the  same  as  the  estimator  obtained  in 
the  Gaussian  analysis  with  equivalent  noise  2q  (a+i^+i.^)  . 
Therefore,  if  F(*)  =  -  0/o^,  the  MAP  estimate  is  after 
Eq.  (4-39) 


p  (r) 


for  i, 
b 


>  >  b  , 


kbo0  k 

^ - :  ~ — ; — r  {-  I  3  .  r  .  cos  (■  ( r ) 

27<(a,-rb+ld)  j  =  1  j  3 


Z  u . r . sinO  (r ) )  , 

] — 1  1  3 


(5-20) 


and  the  ML  estimate  is  the  same  one  'riven  by  I'.q .  (4-48)  . 

In  order  to  illustrate  how  Eqs.  (5-12)  and  (5-17)  change 
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with  k  and  to  compare  them  to  Eqs.  (5-20)  and  (4-48),  some 
examples  are  given  in  Tables  IV  and  V. 


Performance  of  the  Phase  Estimators  in  Shot  Noise 

Considering  that  the  phase  estimates  found  in  the 
Gaussian  noise  analysis  are  asymptotically  unbiased  for 
high  SNR,  the  Cramer-Rao  inequality  seems  to  be  also  an 
appropriate  measure  of  performance  of  the  shot  noise  phase 
estimator  under  the  same  SNR  restrictions.  From  Eqs.  (4-78) 
and  (5-6) ,  the  MAP-CRB  bound  is  given  by 


CRB  = 


-  J, 


k 

I 

Lj=i 


JS  ■-.«» 


1  \-l 


3F(0) 


s .  (0)  + 


{Ab+Xd)T/k 


89 


.  (5-21) 


By  making  the  substitutions  of  Eqs.  (3-25)  and  (4-85)  into 
Eq.  (5-21) ,  the  bound  transforms  into 


CRB  = 


bT 


H 

-  Et 


k 

Z 


[a  jSin9+B  jCos9 ]  2k/2ir 


|  |^27rq  [a  jCOs6-B  j  sin6  ]  k/2u+  (a+i^+i^)  /b  J 

"3F(6)1  l 


30 


I 


(5-22) 


For  Gaussian  phase,  Eq.  (5-22)  becomes 


CRB  =  !  E 


|  2irq 

A  i  bT.  +  JLl' 

\2iq  Ao  2 f 


_ _ _ _  +  1_) 

.  L ,  (a  .cos0-B  ■  sin9)  +2tr  (a+i.+i  .)  /k.  I  2  1 

b=l  3  3  b  d  bj  Oq  ) 


k 

Z 


(g^sinQ+B jCos9 ) 


-1 


9 

1  +  a2bTA  /27rq 
0  o 


(5-23) 


80 


TABLE  IV 


Examples  of  the  MAP  Phase  Estimator  Form 
Under  Multimode  Shot  Noise  Conditions 


DC  Blocked,  Eq.  (5-17),  F(0)  =  -  0/a^ 


D 

A 

6  (i^  ...  rk) 

3 

Cg{r1cot0(r)  -  (r2+r3) tan (0 (r)  +  } 

D 

(rl+r3) oote (r)  -  (r2+r4) tan0 (r) } 

CO 

T 

Og  f  dt  r (t) cot (iot+0) 

Background  Limited,  Eq.  (5-20),  F(0)  =  -  0 /Oq 


D 

®  (^2  •  •  • 

3 

D 

/2  (r1-r3)k1cos§ (r)  -  /2  (r2~r4)k1sin§ (r) 

00 

r  T  ^  -  T  ~ 

[J  r  (t) cosutdt] k  cos0-[J  r  (t) sinutdt] k  sin© 

0  q  0  q 

TABLE  V 


Examples  of  the  ML  Phase  Estimator 
Under  Multimode  Shot  Noise  Conditions 


DC  Blocked,  Eq.  (5-17),  F(0)  =  0 


B 

tan§ (r1  . . .  r^) 

3 

r^  /3  -  tan§ (r) 
r2+r3  /3  +  tan§ (r) 

4 

V  r2  +  r4 

OO 

T 

f  r  (t)  cot  (uit+0  )  dt  =  0 

O 

Background  Limited,  Eq.  (4-48) 


B 

tan§  (r^  .  . .  r^) 

3 

/3  2rl_r2_r3 

3  r2-r3 

4 

rl  ‘  r3 
r2  "  r4 

OO 

J 

J 

rT 

'  r(t)coswtdt 

0 

rT 

r (t) sinwtdt 

0 
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tanQ  (r1  ...  rfc) 


r^/ (sin6 (r) +1. 11)  -  r3/ (sin§ (r ) -1. 11) 


r2/  (cos§ (r) +1 . 11)  -  r^/ (cos6 (r) -1. 11) 


/  dt  r  (t)  cot  +  J)  =  0 


I 


where  Aq  is  the  expected  value  of  the  sum  with  respect  to  9. 
For  (a+i^+i^)  =  b  and  k  very  large.  Eg.  (5-23)  reduces  to 
(Ref  2:300) 


V  (e  )  > 

ar'  r' 


1  +  Oq bT/q 


The  ML  bound  is  on  the  other  hand  (a^  large) 


(5-24) 


Var(er}  -  M  ' 


(5-25) 


the  reciprocal  of  the  collected  count  energy.  The  perform¬ 
ance  of  Eq.  (5-20)  as  given  by  the  CRB  bound  is,  from  Eq. 
(5-22)  with  F ( 0 )  =  -  0/0g, 


CRB  = 


ibT  J 

_  o 

k  (a  .  sin0+6  .  cos0 )  k/2Tr 

y  3  3 

♦  if1 

1 2irq  E 

_j=l  (a+ib+id)/b 

(5-26) 


Invoking  the  identities  of  Eqs.  (4-29)  and  (4-30) ,  Eq. 
(5-26)  simplifies  to 

,-l 


CRB 


.  .< 


b^T  k  ^  2  1  I 

A  +  #  ■ 


(5-27) 


Finally,  the  definitions  of  Eqs.  (4-59),  (4-60)  and  (4-77) 
can  be  used  to  write  the  error  variance  as 


E[e‘]  ;> 


1  +  Eao2/[2TT2q(a+ib+id)Yk] 


(a+ifa+id)  >>  b  , 


(5-28) 
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or  in  the  limit  as  k  becomes  very  large 


l.i.m.  E[eJ]  >  - - -  ,  i  >>  b  .(5-29) 

k  -*■  ®  1  +  E^/q  (a+ib+id)  ^ 


Because  of  the  requirement  that  i^  >>  b,  small  SNR,  the 
error  is  limited  by  the  variance.  Furthermore,  the  actual 
error  variance  may  be  much  higher  due  to  the  inherent  SNR 
restrictions  needed  to  apply  Eq.  (5-22) .  For  the  ML  esti¬ 
mate,  Eq.  (5-28)  becomes 

2 

2tt  q  (a+i,  +i  .) 

Var(er)  - £—2 -  ,  ib  >>  b  ,  (5-30) 

a 

and  Eq.  (5-29)  becomes 


l.i* m . 

k  -*■  °° 


Var<£r> 


q(a+ib+id) 


ib  >  b 


(5-31) 


Equation  (5-31)  tells  that  the  ML  estimate  obtained  from 
Eq.  (4-48) ,  when  used  as  an  approximation  for  background 
limited  conditions,  results  in  a  bad  estimate. 

The  foregoing  discussion  is  based  on  Wyant's  assertion 
(Ref  16:2624)  that  for  shot  noise  limited  conditions,  a 
better  method  for  measuring  phase  is  as  given  by  Eq.  (5-19) , 
F (0 )  =  0,  easily  implemented  by  letting  k  =  4.  (See  Table 
V,  Background  Limited.)  The  performance  of  Eq.  (4-48) 
analogous  to  Eq.  (5-19)  for  ML,  was  proven  to  be  exclusively 
SNR  dependent.  The  equivalent  noise  2q (a+i^+i^)  also  makes 
the  result  of  Eq.  (5-19)  SNR  dependent,  and  the  estimate  is 
likely  to  be  erroneous.  (See  Eq.  (5-31).)  A  better 
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t 


equation  is  given  also  in  Table  V  for  y  -  1,  but  the  solu¬ 
tion  is  not  so  straightforward.  The  analysis  performed  by 
Wyant  gives  an  error  variance  (dark  current  neglected) 


Varur)  -  W)2 


(5-32) 


4y  I  r. 
j  =  l  D 

where  the  photon  counts  can  be  evaluated  from  external 
parameters.  Comparison  of  Eq.  (5-32)  to  Eq.  (5-30)  using 
Eq.  (5-9)  gives  the  following  inequality: 

ir2(a+ib)  2  47r2q(a+i.) 

- k -  -  ~ /  (a+ifa)  >>  b  .  (5-33) 


4b  Z  r. 
j=l  3 


0  9 

d  T  k  Z  a? 
3=1  3 


For  the  particular  case  of  k  =  4,  Eq.  (5-33)  yields 


(a+ib)  ;>  |(r1+r2+r3+r4)  >>  b  . 


(5-34) 


From  Eqs .  (5-32)  and  (5-9),  Wyant’s  approach  is  found  to  be 

also  SNR  dependent.  Any  good  performance  thus  depends  on 
the  condition  that 


£  r.  >>  (a+i  ) 
j=l  3  b 


(5-35) 


From  Eq.  (5-34),  obviously  Eq.  (5-35)  is  not  true. 

The  foregoing  results  clearly  indicate  that  phase  esti¬ 
mation  requires  a  strong  signal  regardless  of  the  noise 
process.  It  can  be  noted  that  under  low  light  level  condi¬ 
tions,  the  SNR  constraints  are  more  difficult  to  meet,  and 
large  errors  may  be  expected.  With  these  conclusions,  the 
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analysis  of  phase  difference  estimation  from  measurements 
from  a  single  detector  are  completed  within  the  scope  of 
this  thesis.  Given  the  added  complexity  of  Eq.  (5-8)  and 
time  constraints,  a  simulation  as  performed  in  Chapter  IV 
is  not  included.  Estimation  of  the  actual  phase  components 
in  Eq.  (2-3)  will  be  performed  in  the  next  chapter  using 
the  measurements  of  the  detector  array  as  a  whole. 
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VI  Joint  Processing  of  Array  Signals 
for  Wavefront  Estimation 

The  problem  of  phase  difference  estimation  was 
addressed  in  Chapters  IV  and  V.  Ideally  what  is  desired 
are  the  actual  phases  at  specific  points  across  the  aper¬ 
ture  rather  than  slopes  between  points.  The  methods  usually 
employed  to  find  these  phases  consist  of  mappings  of  data 
given  by  the  measured  wavefront  difference  functions  A<J>  (r  ) 
and  the  geometry  of  the  data  points.  These  mappings  use  the 
concepts  of  least  squares  fitting  and  are  discussed  by 
Fried  (Ref  1) ,  Hardy  et  al  (Ref  4) ,  Hudgin  (Ref  5) ,  Rimmer 
(Ref  10),  and  Wyant  (Ref  16)  among  others.  The  general 
result  given  by  these  mappings  is  an  average  of  phases  and 
measurements  about  a  single  point,  requiring  a  recursive 
solution. 

A  different  approach  will  be  used  in  this  thesis  which 
uses  the  fact  that  the  phases  across  the  aperture  are  spa¬ 
tially  correlated  due  to  slow  spatial  variation  of  the 
wavefront,  and  the  assumption  that  a  spatial  covariance 
matrix  is  available  from  experimental  measurements.  There¬ 
fore,  by  jointly  processing  the  outputs  of  the  two  detector 
arrays  discussed  in  Chapter  II,  a  phase  estimate  can  be 
made  in  real  time.  Because  of  the  a  priori  information 
supplied  by  the  covariance  matrix,  an  improvement  is  expec¬ 
ted  over  mappings  of  data  points. 

The  criteria  of  Maximum  A  Posteriori  estimation  will 
be  used  in  this  chapter,  but  the  algorithm  so  obtained  will 
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be  in  continuous  waveform  rather  than  discrete  counts. 
Because  of  the  difficulty  in  applying  the  MAP  theory,  sev¬ 
eral  assumptions  need  to  be  made  in  order  to  simplify  the 
algorithm  derivation.  The  problem  will  be  worked  out  in 
the  Gaussian  context  or  second  moment  models  only,  where 
detector  thermal  noise  is  the  predominant  noise  source.  The 
noise  waveforms  from  each  detector  are  samples  from  indepen¬ 
dent  white  Gaussian  processes  with  zero-mean  and  strength 
N^/2.  The  subscript  k  will  be  used  to  index  the 
detector.  The  phases  <j>(r,t)  will  be  assumed  to  be  stepwise 
constant  in  the  interval  sequence  (0,  T,  2T,  .  ..)  as  pre¬ 
sented  in  Chapter  I.  The  covariance  matrix  is,  therefore, 
constant  in  each  interval ,  but  needs  to  be  updated  every  T 
seconds.  This  is  the  sequential  problem  which  is  beyond 
the  scope  of  this  thesis.  A  suggestion,  however,  will  be 
given  later  for  sequential  estimation  by  quantizing  0(r,t) 
when  it  is  continuous  in  time.  The  final  assumption  to  be 
made  is  that  of  a  jointly  Gaussian  random  phase  distribution 
over  the  aperture.  The  time-space  problem  fitted  to  the 
above  description  is  one  of  multiple  channel,  multiple 
parameter  estimation  in  Gaussian  noise. 

Multidimensional  Estimator  Formulation 

Derivation  of  the  required  time-space  estimator  equa¬ 
tion  will  be  made  in  this  section.  The  next  section  will 
treat  the  specific  application  to  the  shearing  interferome¬ 
ter  output.  So,  the  multiple  observation  model  can  be 
written  in  vector  form  as 


r(t)  =  s(t,a)+n(t),  0  <  t  <  T  , 


(6-1) 


where  r(t)  is  a  column  vector  of  the  outputs  from  a  two- 
dimensional  detector  array,  indexed  with  a  single  subscript 
1  <  k  <  m.  A  vector  of  phase  parameters  a  over  the  aper¬ 
ture,  indexed  also  with  a  single  subscript  1  <  p  £  n,  is  to 
be  estimated  jointly  using  all  available  outputs.  The  for¬ 
mulation  needed  for  array  processing  is  a  direct  extension 
of  the  single  element  case.  Let  the  output  from  the 
detector  be 


r  K(t)  =  sK(t,a)  +  nK(t) 


(6-2) 


The  noise  statistics  are  given,  from  Eqs.  (4-2)-(4-6),  as 


E[nK(t)]  =  0 


E[nK(t)nK(t')]  =  -ir^  S(t-t') 


(6-3) 


*  W 
A  K 


6 (t-t ' )  . 


(6-4) 


It  is  assumed  that  the  noises  in  the  detectors  are  statis¬ 
tically  independent  of  each  other.  The  detector  outputs  r 
are  also  statistically  independent  as  was  shown  in  Eqs. 
(4-9) - (4-13) .  Therefore,  an  array  output  vector  can  be 
defined  as 


R  -  [£x  . 


.  r  . .  .  r  ] 

-<  -m 


(6-5) 


i 


and  the  conditional  density  of  the  observations  given  the 
phase  vector  is  found  from 


m 

n 

K  =  1 


(r  la) 


-K  1  — 


(6-6) 


where  f  ,  (r  |a)  is  given  by  Eq.  (4-14'  The  likelihood 
— k  *  — 

function  defined  by  the  ratio 


A 


A 


(6-7) 


is  therefore  given  by 


A  =  ^k 
T 


m 

Z 


k 

Z 


K=1  j=l 


r  . 
1 


(>  ) 


W 


s\<]  (*}  "  2W~  (sjK)  (2:))2 


(6-8) 

Equation  (6-8)  must  be  maximized  by  the  proper  choice  of 
all  the  elements  a^  in  a.  Since  the  phases  across  the  aper¬ 
ture  are  assumed  jointly  Gaussian  and  spatially  correlated, 
they  can  be  represented  in  a  different  coordinate  system 
where  the  new  elements  a^  are  independent  Gaussian  random 
variables,  each  with  density 


f  (a.  ) 
a .  i 
l 


[2tt o2  ]  1//2exp  [-a?/2 a2  } 

d  .  -L  ct  . 

i  i 


(6-9) 


Beginning  the  estimator  derivation  with  Eq.  (6-9) ,  the  MAP 
estimate  of  a.  can  be  obtained  by  maximizing  {A  +  Inf  (a.)} 

1  cl  .  1 

1 

with  respect  to  a^.  The  result  in  vector  form  is 


_  2,  k  8s1  (a) 

2q  jj  y  3 
T  3a^  -  l-j  -j 


W  1 [r .  -  s. (a) ]  - 


a . 

l 


=  0 


(6-10) 
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By  using  Eqs.  (6-14)  and  (6-15),  Eq.  (6-10)  becomes 


a .  = 

l 


o2  I  <J)TD,  (a)W_1[r .  -  s.(a)] 
T  ai  j=i  -i-D  -  ~3  ~3  ~ 


(6-16) 


By  using  Eq.  (6-16)  in  Eq.  (6-12)  ,  the  joint  estimate  is 


a  = 


2q2k 


N  2  T 
1 .  l .  m.  Z  a  <b  .  dr. 

N  -*■  00  i=l  ai  -1-1 


,-l 


l  D. (a) W  [r  .  -  s .  (a) ] 
j  =  l  “3  _  “3  ~3  ~ 


The  covariance  matrix  is  given  by  (Ref  12:222) 

N  2  T 

k  =  l.i.m.  Z  cr  <f> .  <t> . 

N  -*■  00  i=l  ai  1  1 


(6-17) 


(6-18) 


and  is  assumed  to  be  known  in  the  interval  (t^t^T).  It 
is  further  defined  by 


k  -  E [a  aT] 

cl 


k,  .  .  .  k, 

1  In 


k  .  k 

nl  n 


(6-19) 


Therefore,  the  discrete  joint  MAP  estimate  of  a  is  given  by 

a  (r)  =  k  Z  D .  (a)  W  [r  .-s  .  (a)  ]  ,  t.  £  t  <  t.+T  . 

—  —  x  — a  —  j  —  —  —  j  ~J  —  J- 

(6-20) 

Equation  (6-20)  can  be  converted  into  a  continuous  form  by 
substituting  for  the  definition  of  Eq.  (3-14) .  By  observing 
that,  for  a  set  £ (t)  of  complete  orthonormal  functions, 
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a 


6 (t-u) 


(6-21) 


k 

l.i.m.  Z  5. (t)c-(u)  = 

k  ->•  00  j  =  l  “*  -* 

then  Eq.  (6-20)  becomes 


a  =  2k 


ka  J  D (t , a) W-1 [r 


(t)  -  s (t ,a) ] dt 


for  0  <  t  <  T  , 
where 


(6-22) 


D  (t ,  a)  - 


Zs1 (t,a) 
9a, 


(t,a) 
3a 


9sm(t,a) 

dal 


9sm(t,a) 


(6-23) 


3a 


n 


Equation  (6-22)  is  the  joint  MAP  estimate  from  a  continuous 
measurement  of  a  Gaussian  random  vector  a  observed  in  Gauss¬ 
ian  noise.  This  is  the  result  that  will  be  used  in  the 
forthcoming  developments  to  process  jointly  the  outputs  of 
a  shearing  interferometer.  Equation  (6-22)  can  also  be 
obtained  in  a  manner  analogous  to  the  single  detector  by 
performing  a  correlation-summiny  operation  as  shown  in 
Figure  20  (Ref  14:367,452,453). 


Wavefront  Estimation  from  the  Shearing  Interferometers 

The  output  fields  of  the  interferometers  at  the  focal 
planes  are  received  and  processed  by  two  separate  detector 
arrays,  one  each  for  the  X  and  Y  sheared  fields  as  depicted 
in  Figure  2.  Therefore,  it  will  be  convenient  to  keep 
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Fig  20.  Correlator-Summer  Model  for  Joint  Processing 

track  of  the  observations  by  use  of  appropriate  subscripted 
notation.  Define 


tor  in  each  of  the  x-  and  y-arrays.  The  subscripts  <  and  l 
index  the  rows  and  columns  respectively.  For  any  one  detec¬ 
tor,  the  signals  can  be  written  as 

y(t,Q)  =  a  +  b  sin  (wt  +  8  (y  ) ) 


for  fixed, 

a 


(6-26) 


l 


and 


x(t,0)  =  d  +  c  sin  (tot  +  6  (x  )) 

cl 


for  v  fixed, 

ci 


(6-27) 


where  the  difference  functions  are,  from  Eqs.  (2-3)  and  (3-1) 
(3-5) ,  given  by 

9  (x  )  =  \  (4>(x  -  M  ,y  )  -  <j>(x_  +  M  ,y_)}  (6-28) 


a  sd  a 


a  s, '•*  a' 
a 


0  (ya)  =  \  ^ (xa » ya  -  Ms  )  -  <}>(xa,y  +  Mg  )}  .  (6-29) 

d  a 

Further  notational  simplification  can  be  made  by  defining 


<J>  (x“) 

A 

*(xa  '  Msd'ya> 

(6-30) 

4>  (x+) 

A 

<Mxa  +  MSd,ya) 

(6-31) 

(y~) 

A 

*(xa'ya  ~ 

(6-32) 

+~~ 

-e- 

A 

^  ^xa,ya  +  MSj}  ' 

(6-33) 

so  that  Eqs.  (6-26)  and  (6-27)  can  be  written  as 


y(t,0)  =  a  +  b  sin(ut  +  — 2  —  - ) 


for  x  fixed, 
a 

and  x(t,0)  =  d  +  c  sin  (cot  +  -1 — —  ■  ) 


for  y  fixed, 
a 


(6-34) 


(6-35) 


Equation  (6-22)  can  now  be  applied  using  Eqs.  (6-24)  through 
(6-35)  to  perform  the  joint  estimate  of  the  wavefront  phases 
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using  two  plane  detector  arrays  properly  interfaced  to 
couple  the  x-  and  y-measurements  with  the  same  wavefront 
points.  But  the  two-array  configuration  must  be  set  up 
before  proceeding  to  perform  the  joint  phase  estimation. 

Configuration  of  the  Detector  Arrays.  The  need  for  the 
use  of  two  detector  plane  arrays  for  actual  phase  estimation 
is  a  consequence  of  the  structure  of  the  difference  func¬ 
tions  in  Eqs .  (6-30)  through  (6-35),  obtained  from  the  use 
of  lateral  shearing  instead  of  radial  shearing  interferom¬ 
etry.  As  is  the  case  in  the  approach  of  data  mappings,  pro¬ 
cessing  of  the  observations  r(t)  requires  a  specific  detec¬ 
tor  arrangement.  In  order  to  make  full  utilization  of  the 
information  collected,  the  detectors  (and  the  arrays)  must 
be  arranged  so  that  each  wavefront  point  be  measured  by  as 
many  detectors  as  possible  in  order  to  provide  a  strong 
deterministic  relationship  between  the  m >asurements  over  the 
aperture.  All  other  coupling  is  provided  statistically  by  the 
covariance  matrix  of  the  phases. 

In  order  to  provide  redundancy  of  measurements,  the 
best  possible  detector  arrangement  is  as  shown  in  Figure  21. 
This  arrangement  allows  each  phase  point  to  be  measured  by 
four  detectors,  two  from  each  one  of  the  arrays,  and  has 
the  advantage  that  only  one  reference  phase  is  required  to 
determine  the  entire  phase  distribution.  It  is  assumed  that 
such  phase  point  is  measured  by  another  means  or  is  set 
arbitrarily  equal  to  zero.  The  X's  and  Y's  in  Figure  21 
denote  the  detector  locations  and  the  4>'s  denote  the  phase 
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Fig  21.  Arrangement  of  Two  Detector  Arrays  to  Measure 
(m+l)x(n+l)  Phase  Points  with  (m+l)xn  X-  and 
mx(n+l)  Y-Detectors 


points  being  measured.  Except  for  the  edge  phase  points, 
each  phase  <PK  ^  is  measured  by  detectors  xk£,  Y  k%' 

and  y^  „  .  ,  4>  being  a  common  phase  point  of  the  four  dif- 
ference  functions  0K*(xa),  0K_1# A  (xa> ,  e<£(ya>,  and  6(c^_1(ya) 
The  grid  shown  in  Figure  21  is  formed  by  overlapping  the  two 
arrays  to  indicate  the  relative  positions  of  the  X-  and  Y- 
shear  detector.  The  corresponding  detectors  x^  and  y^ ^  are 
not  located  at  the  same  point  on  the  field.  They  are  dis¬ 
placed  45°  instead  so  that  both  can  measure  the  same  phase 
^ .  The  direction  of  the  shear  has  been  chosen  from  right 
to  left  and  bottom  up  to  correspond  to  the  notation  adopted 
in  Eqs.  (5-34)  through  (6-37).  In  this  arrangement  the 
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detector  spacing  has  been  chosen  so  that 


=  *K+l,£(y+) 

(6-36) 

and 

*k£(x_) 

=  *k,£+1{x+) 

(6-37) 

in  Eqs.  (6-30)  through  (6-35) .  The  detectors  must  then  be 
spaced  2M  shear  distances  s^  apart  (s^x  and  s^  may  not 
necessarily  be  equal) .  This  is  depicted  in  Figure  22  for  a 
column  array. 

In  general,  for  a  non-square  array  of  phase  points,  the 
number  of  detectors  is  different  in  each  array.  Thus, 
referring  to  Figure  21,  there  are  [m+l]xn  X-detectors  and 
mx[n+l]  Y-detectors,  and  [m+1] [n+1]  measured  phases.  The 
ratio  of  detector  to  phase  points  is 


[m+1] xn  +  nx [n+1] 
[m+1] x [n+1] 


and  ranges  from  unity,  when  a  minimum  of  four  phases  are 
measured  with  two  each  X-  and  Y-detectors,  up  to  the  value 
of  two,  in  the  limit  when  a  very  large  number  of  phase 
points  are  being  measured.  Thus,  in  the  best  possible 
arrangement  of  Figure  21,  the  number  of  detectors  required 
tends  to  double  as  more  wavefront  points  are  measured  simul¬ 
taneously.  With  these  preliminaries  completed,  derivation 
of  the  algorithm  for  wavefront  estimation  by  jointly  pro¬ 
cessing  the  two  detector  array  outputs  can  be  initiated. 

Joint  Wavefront  Estimation  with  Two  Detector  Arrays. 

The  observations  r(t)  can  be  arranged  in  a  column  vector  of 
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Fig  22.  Arrangement  of  Detector  Column  Array  for 
Joint  Processing 


dimension  (m[2n+l]+n),  using  the  notation  of  Eqs.  (6-24)  and 
(6-25) .  Thus , 


— (t)  [zll (t) 'wll (t)  zmn(t) ,wmn(t) 


•  Wm+l,n(t)] 


(6-3f 


where  z^ft)  and  wK^(t)  are  the  observations  with  signal 
components  given,  after  Eqs.  (6-34)  through  (6-37),  by 

4>. 

for  0  <  t  <  T, 

1  s  k  <  m, 

1  <  S,  <  n+1. 


Y^o(t.i)  =  +  bKAsin(wt  +  — +1'-^ - — ] 
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(6-39) 


£+1 

and  =  d<£  +  sin(wt  +  — 1 - ^ ~ > 

for  0  <  t  <  T, 

1  <  k  <  m+1, 

1  <  £  <  n  .  (6-40) 

The  arguments  (x+)  and  (y+)  used  in  Egs.  (6-34)  and  (6-35) 
are  dropped  from  Eqs.  (6-39)  and  (6-40)  since  the  ambiguity 
is  taken  care  of  by  the  subscripts  k+1,£  and  k,£+1,  and 
because  <(>K)i(x+)  =  <J>K^(y+). 

The  phase  vector  is  also  a  ([m+1] [n+1] )  column  matrix 


=  [<J> 


11 


pq 


^m+1 ,n+l 


(6-41) 


and  has  a  ([m+1] [n+1])  symmetrical  covariance  matrix  with 

terms  from  kU;U  to  Vl ,„+l ;m+l  ,n+l •  The  noise  ^variance 
matrix,  on  the  other  hand,  is  (m[2n+l]+n)  diagonal  (noises 
were  assumed  spatially  uncorrelated)  given  by 


W 


M 


11 


O 


O 


n,  n+1 


(6-42) 


|_  Mm+l,nJ 

where  the  N's  and  M's  make  reference  to  X-  and  Y-detector 
noises  respectively.  With  the  observation  and  covariance 
matrices  defined  in  Eqs.  (6-38)  through  (6-42) ,  Eqs.  (6-22) 
and  (6-23)  can  be  applied  directly  to  perform  the  joint 
estimation  of  the  vector  <p. 
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Thus,  the  product  W  ‘*'[r(t)  -  s(t,<}0]  in  Eg.  (6-22)  is 
a  (n[2n+l]+n)  column  matrix  with  general  terms 

N~  t ZKZ(t)  ~ 

K  Jc 

for  1  <  <  <  m, 

1  <  £  <  n+1,  (6-43) 

and 

for  1  <  k  <  m+1, 

1  £  Z  <  n  .  (6-44) 

The  signal  derivative  matrix  is  ( [m+1] [n+1] ) x (m [ 2n+l] +n) 
dimensional  given  by 


the  MAP  estimate  of  the  pq 
terms  to  be  given  by 


phase  is  found  in  general 


3y^(t,i) 


m+l 

Z 

i=l 


m  n+1  2k 
Z  Z  — 
k=1  £=1 


EZiil 

N.,n 


3<J> 


ID 


[ZKZ (t) ~YKl (t'i} ]dt 


m+l  n  /*T  3x  ,  (t,40  ) 

+  E  I  /  - - [w  (t)-x  (t,<!»)]dt  .  (6-46) 

k=1  Jl=l^  d$ij  "  ) 

Finally,  when  Eqs.  (6-39)  and  (6-40)  are  used  in  Eq.  (6-46) , 
the  general  pq^^  member  of  the  <)>_  vector  of  Eq.  (6-41)  is 
obtained.  The  result  is  given  by 


m  n+1  b 


^Pq  Kh  NkJI  tkpq,*K+l,£  "  kpq,K£] 

T  ~ 

.  ,  .  ,  ^K+l,£ 

zKji  (t)cos  (cot  +  - '-~2 - )  dt 


m+l  n  C 


kZ 


+  ^  ^  m  I^CCI.  k:  9,+l  ~  ^na 

k=1  £=1  k£  Pq'K^+1  Pq,KJt 


J  WK  (t)cos(wt  +  \ - — )dt 


for  0  <  t  £  T, 

1  <  p  <  m+l , 

1  £  q  <  n+1,  (6-47) 


Equation  (6-47)  is  the  joint  Maximum  A  Posteriori  phase 
estimator  using  measurements  of  two  orthogonal,  lateral 
shearing  interferometers  and  their  detector  arrays.  The 
phase  estimate  distribution  over  the  aperture  is  shown  in 
Figure  23. 


1.  J 


Fig  23.  Phase  Estimate  Distribution  Over  the  Aperture 

The  algorithm  given  by  Eq.  (6-47)  is  applicable  only 
when  the  phases  <£  being  measured  are  jointly  Gaussian  ran¬ 
dom  variables  and  the  a  priori  information  represented  by 

k_  is  available.  The  distinction  between  <J>  and  $  in 
—a  pq  rrs 

Eq.  (6-47)  is  made,  observing  the  structural  form  of  Eq. 

(6-46),  by  the  covariance  terms  k  .  .  and  k  .  .  only, 

pq  r  ij  rs  ,  ij 

where  i  and  j  take  on  all  values  from  1  to  m+1,  and  1  to 
n+1  respectively.  Therefore,  the  joint  phase  estimates  are 
weighted  accordingly  by  the  a  priori  information.  This 
will  be  further  explained  in  a  forthcoming  example. 

Solution  of  Eq.  (6-47)  is  to  be  obtained  recursively 
with  numerical  methods  on  a  digital  computer.  However,  to 
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illustrate  the  estimator,  Eq.  (6-47)  can  be  implemented  with 
a  heuristic  correlator-summer  of  considerable  complexity  as 
will  be  seen  in  the  following  example. 


Illustrative  Example.  The  simplest  example  to  illus¬ 
trate  a  hardware  implementation  of  Eq.  (6-47)  is  given  by 
the  joint  estimate  of  a  wavefront  at  four  locations  using 
two  (2x2)  detector  arrays.  The  arrangement  is  shown  in 
Figure  24.  The  notation  is  so  chosen  for  simplicity. 


Detector  Y2  corresponds  to  and  phase  corresponds  to 

<t>09  in  the  notation  of  Eq.  (6-47).  It  will  also  be  assumed 


that  the  noises  are  of  equal  strength  Nq  and  the  signals 
have  equal  amplitude  /2Ea/T.  With  the  problem  so  defined, 

/  i  u  \ 

the  p'  phase  estimate  is  given,  from  Eq.  (6-47),  by 


0 

P 


^  ( 


f  4*  i 

(kp3-kpg)  /  Zj^CtJcosfut  +  — ~^)  at 


0Y  T  p3  pi 


r  v*2 

+  (k  4~k  2)  /  (t)  cos  (cut  +  — ^ — ~)  dt 


T 

w  (t)  cos  (cot  +  - — 2 — dt 


r  V'N  I 

+  (k  4~k  3)  /  w2(t)cos(ot  +  -  2~  -)  dt .  ,  (6-48) 

J0 


where  any  <t>  is  set  as  the  zero  reference.  Equation  (6-48) 
can  be  implemented  with  a  correlator-summer  as  shown  in 

O' 

4> .  -4>  . 

Figure  25,  where  cos  (wt  +  1  ^ )  is  denoted  by  ^  The 

z  r  1 

structure  of  Figure  25  is  a  feedback  system  in  which  the 
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Fig  24. 

Phase-Detector 

Grid  for  Joint 

Processing  of  Four 

Detectors  to  Estimate  Four  Wavefront  Phases 


detector  outputs  are  weighted  by  the  covariance  terms  so 
that  each  phase  estimate  has  a  contribution  from  each  meas¬ 
urement.  Three  feedback  loops  from  the  estimate  outputs 

A  A  A 

<t>2,  $3  and  <J>4  are  returned  to  the  detector  inputs.  The 
estimate  <J>^  does  not  provide  any  feedback  for  being  the 
starting  point  of  the  spatially  recursive  estimation.  The 
phase  output  $3  has  been  arbitrarily  set  as  the  zero  refer¬ 
ence. 

To  better  understand  how  this  formidable  structure 
weighs  the  measurements,  or  better  yet,  what  Eq.  (6-48) 
really  does,  the  underlying  mechanism  is  shown  in  Figure  26. 
Here,  the  phase  points  have  been  replaced  by  the  covariance 
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Fig  25.  Correlator-Summer  to  Estimate  Four  Phases 


I 


terms  of  <j>  with  respect  to  the  other  phases.  The  detector 
locations  have  been  denoted  by  their  measurements  z(t)  and 
w(t).  Figure  26-a  shows  the  mechanism  involved  in  esti¬ 
mating  Each  phase  location  has  the  covariance  term 

involving  <J>^  and  the  phase  at  each  point;  this  is  the  mean¬ 
ing  of  the  arrows.  The  operation  involves  taking  pairwise 
the  difference  between  the  covariance  terms  at  all  points 
and  multiplying  by  the  measurement  made  between  each  two 
points  (this  means  the  output  of  the  integrator) .  Thus, 
the  difference  between  and  k^4  is  multiplied  by  zero 

because  there  is  no  detector  in  between.  In  the  same  man¬ 


ner,  Figure  26-b  shows  the  same  mechanism  used  to  evaluate 
4>3.  Here,  the  covariance  terms  are  k33,  k31,  ^32  and  k34* 
This  can  be  generalized  for  any  number  of  array  elements. 

From  this  example,  it  is  evident  that  the  algorithm 
indeed  provides  a  means  for  jointly  processing  spatially 
correlated  phase  measurements.  It  also  implies  a  simultane¬ 
ous  evaluation  of  all  phase  points  and  a  recursive  substi¬ 
tution.  The  difficulty  involved  in  solving  Eq.  (6-47) 
numerically  is  also  evident  due  to  the  redundancy  provided 
by  the  feedback  loops.  Therefore,  such  a  solution  goes 
beyond  the  scope  of  this  paper. 


Estimator  Performance 

The  performance  of  the  estimator  of  Eq.  (6-22)  is  more 


easily  described  in  terms  of  the  Cramer-Rao  bounds,  provided 
that  the  signal-to-noise  ratio  is  large  and  biases  are 


negligible.  Then,  the  lower  bounds  on  the  error  matrix  are 
the  diagonal  terms  of  the  matrix 

SB  -  5a(I  +  W1  (6-49) 

adapted  from  a  more  general  case  given  by  Van  Trees  (Ref  13: 
454).  Simplification  to  the  point  of  Eq.  (6-49)  is  possible 
because  the  matrix 

D_(u,4>)  - 

— a  — 

A 


T 

/ 


E[D(t,<f.)W"1DT(t,40  ]dt 


°a 


(6-50) 


is  constant  over  (0,T)  for  s(t,£)  defined  as  [a+bsin (wt+£) ] . 

To  find  the  MS  error  bounds,  the  matrix  operations  of 
Eqs.  (6-50)  and  (6-49)  must  be  performed.  Carrying  out 
these  matrix  multiplications  is  a  rather  cumbersome  task. 
After  some  work,  the  matrix  D  (u,<|>)  is  found  to  have  a 

—a  — 

diagonal  band  form  with  dimensions  (m[2n+l] +n) x ( [m+1] [n+1] ) 
as  illustrated  in  Figure  27.  The  five  X's  in  each  column 
are  niven  by 


2  /3yp--V>q  3yp— .Us) 

p-l ,q  Wp-l,q  3°pq  / 

*p,  q-1  Wp,q-1  3*pq  / 


(6-51) 


(6-52) 
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The  product  D  (u,<+>)k  is  a  ([m+l][n+l]  )  square  matrix 

—cl  —  — a 

general  term  at  location  pq,ij  is  given  by 


< 


for 

N,y:  1  <  P  <  m  ,  1  <  q  <  n+1 

M,x:  1  <  p  <  m+1  ,  1  ^  q  <  n 

1  <  i  i  m+1  ,  1  $  j  <  n+1  .  (6-56) 

Integration  of  the  matrix  with  general  term  given  by  Eq. 

(6-56)  gives  a  constant  matrix  whose  general  term  at 

location  pq,ij  is  given  by 


2 

-  W  k  Cp,q-1 

4Vl,q  “Mp.g.! 

+  J  +  J~1,<3  +  ia  +  CP.q- 

T  lNpq  Vu  %  «p.o- 


•  M  P-l,q  Mpq  Mp,q-l)  M.ll 

c2  h2 

-  k  pq  . 

Mpq  P '<3+1; ij  ?Npq  Xp+l,q;ij  * 


(6-57) 
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The  matrix  D  is  ([m+1] [n+1] )  square  and  p,  q,  i,  and  j  are 

— Cl 

bounded  as  in  Eq.  (6-56) .  Further  evaluation  of  the  error 
matrix  is  to  be  done  numerically  to  yield  the  diagonal 
values  of  Rg  defined  by  Eq.  (6-49) .  This  evaluation  as 
well  as  a  numerical  solution  of  Eq.  (6-47)  are  not  included 
in  this  paper  due  to  the  difficulties  encountered  in  under¬ 
taking  that  task.  This,  however,  does  not  discourage  fur¬ 
ther  study  since  Eqs.  (6-47)  and  (6-57)  are  believed  to  be 


correct. 


VII  Conclusions  and  Recommendations 


Conclusions 

The  purpose  of  this  work  was  twofold:  (1)  perform  a 
discrete  phase  estimation  analysis  for  a  single  detector 
measurement,  and  (2)  perform  a  joint  phase  estimation 
analysis  for  multiple  detector  measurements.  This  work  was 
prompted  by  the  need  for  improved  processing  techniques 
directly  applicable  to  shearing  interferometry  and  wave- 
front  correction  systems. 

'  The  first  workfront  motivation  was  to  determine  in  a 
stochastic  sense  if  phase  estimation  algorithms  with  the 
simplicity  of  the  form  of 

*  =  tan"1  [§=§]  ,  (7-1) 

from  Reference  16,  intended  for  use  with  low  level  signals, 
could  be  obtained  using  Maximum  A  Posteriori  and  Maximum 
Likelihood  estimation  theories,  and  the  conditions  for  which 
they  would  give  good  phase  estimates.  A  family  of  such 
algorithms  was  found  and  is  given  by  the  ML  estimator  Eq. 
(4-48)  ,  which  is  a  specialized  result  of  the  MAP  Eq.  (4-36) , 
derived  under  the  white  Gaussian  noise  assumption.  A  simi¬ 
lar  result  is  also  given  by  Eq.  (5-8),  derived  under  the 
Poisson  shot  noise  assumption,  for  the  case  of  low  fringe 
visibility.  It  is  shown,  however,  without  empirical  veri¬ 
fication  (computer  simulation) ,  that  this  algorithm  form 
will  result  in  poor  estimates  under  those  modeling  conditions. 
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The  performance  analyses  carried  out  in  Chapters  IV 
and  V  (both  theory  and  simulation)  show  that  the  only  factor 
affecting  phase  estimation  performance  is  the  SNR  regardless 
of  the  noise  process  involved  and  the  technique  used.  In 
fact,  the  performance  of  both  MAP  and  ML  estimators  is 
asymptotically  given  by  the  reciprocal  SNR  with  an  appro¬ 
priate  scaling  factor  for  the  particular  algorithm  form 
used  (given  k) .  Emphasis  is  made  on  this  point  to  clear  up 
Wyant's  implication  that  Eq.  (7-1)  might  be  free  from  the 
SNR  restrictions. 

The  estimators  of  Eqs.  (4-39)  and  (4-48)  were  verified 
with  a  simple  computer  simulation  with  results  presented 
in  Figures  10-19.  The  similarity  of  performance  between 
the  MAP  and  ML  estimators  particularly  for  SNR  >  10  dB 
follows  directly  from  the  implications  of  Eqs.  (4-91)  and 
(4-95) ,  which  predict  the  same  asymptotic  performance  of  the 
MAP  and  ML  phase  estimators  for  large  SNR,  and  the  implica¬ 
tion  of  Eq.  (4-90) ,  which  tends  to  ignore  a  priori  informa¬ 
tion  as  the  noise  in  the  measurement  decreases. 

For  a  given  SNR,  there  exists  a  tradeoff  between  algo¬ 
rithm  simplicity  and  algorithm  performance  of  the  discrete 
estimators.  The  increased  structural  complexity  as  k 
increases  from  k  =  3  to  infinity  is  illustrated  with  a  few 
examples  in  Tables  I,  II,  IV  and  V.  This  complexity  is 
particularly  noticeable  in  the  Poisson  analysis  equations. 
The  return  of  using  the  more  complex  forms  is  an  improved 
performance  as  shown  in  Table  III  and  plotted  on  a  relative 
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basis  in  Figure  9.  The  difference  in  theoretical  perform¬ 
ance  between  the  worst  (k  =  3)  and  the  best  (k  -  »)  possible 
estimator  forms  is  1.65  dB  according  to  the  definition  of 
Eq.  (4-76) ,  and  decreases  quite  rapidly  as  the  number  of 
counts  (k)  increases.  On  the  basis  of  this  performance- 
complexity  tradeoff,  the  photon  counting  technique  (k  <  <*>) 
is  suboptimal  (but  acceptable)  with  respect  to  current- 
measurement  (k  =  co)  based  phase  estimation. 

The  second  motivation  of  this  work  was  the  analysis  of 
a  time-space  problem  intended  to  provide  a  joint  estimate 
of  the  phases  across  the  aperture  of  the  interferometer. 

Such  joint  processing  had  not  been  addressed  in  the  light 
of  MAP  theory,  where  the  fact  that  the  phases  were  spatially 
correlated  could  be  used  to  improve  performance.  The  algo¬ 
rithm  derived  is  given  by  Eq.  (6-47)  and  is  restricted  to 
the  assumption  of  a  Gaussian  phase  distribution.  The 
weighting  between  measurements  is  explicitly  shown  by  the 
covariance  terms  in  Eq.  (6-47) .  The  difficulty  with  this 
algorithm  is  the  mathematical  form  of  simultaneous  non¬ 
linear  integral  equations  for  which  a  solution  is  not 
readily  available.  The  performance  equations  were  carried 
out  up  to  the  point  of  numerical  evaluation,  which  is  not 
included  due  to  time  constraints  imposed  by  the  difficulties 
encountered  in  deriving  the  algorithms.  Thus,  a  direct 
comparison  to  other  estimators  such  as  the  one  derived  by 
Hudgin  (Ref  5)  is  not  possible  given  the  form  of  Eqs.  (6-47) 
and  (6-57). 
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Recommendations 

The  study  contained  in  Chapters  IV  through  VI  was 
performed  for  the  extreme  cases  of  detector  limited  and 
signal  limited  noise  statistics.  It  is  not  unreasonable  to 
think  that  the  intermediate  cases  are  also  encountered  in 
the  reception  of  optical  fields,  where  both  thermal  and 
shot  noise  occur  together.  These  noise  processes  are  inde¬ 
pendent  of  each  other,  and  the  density  of  the  observables 
is  then  a  convolution  of  Gaussian  and  Poisson  functions. 
Estimation  under  these  conditions  should  be  tried. 

The  random  phase  in  the  argument  of  the  sine  function 
represented  both  target  and  turbulence.  Going  beyond  the 
application  of  wavefront  correction  systems,  it  may  be 
desirable  to  distinguish  target  and  noise  induced  phases. 

This  is  estimation  in  the  presence  of  unwanted  parameters 
and  should  also  be  considered. 

The  basic  assumption  of  the  analyses  presented  in  this 
paper  was  the  time  invariance  of  the  phase  in  the  measurement 
interval.  Although  the  staircase  approximation  to  the  phase 
process  may  be  suitable  for  slowly  varying  fields,  it 
requires  an  update  of  the  covariance  matrix  in  each  inter¬ 
val.  A  procedure  should  be  tried  for  which  0 (t)  is  time- 
dependent  throughout  the  measurement,  thus  freeing  the  esti¬ 
mation  problem  from  the  requirement  of  short  observation 
intervals.  This  could  be  done  by  homodyning  the  detector 
signal  and  filtering  the  frequency  domain  components  with  a 
low  pass  filter.  The  output  signal  would  then  have  the  form 
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r(t)  =  B  ~  sinu(t)  +  —  n  (t)  (7-2) 

2  /2  C 

where  B  is  a  gain  factor  given  by  the  amplitude  of  the  mix¬ 
ing  signal,  and  nc(t)  is  filtered  white  Gaussian  noise. 

This  waveform  can  now  be  used  in  Eq.  (3-15)  to  obtain 

r.  =  -=-  smO  .  +  —  n.  .  (7-3) 

3  2  3/J3 

The  signal  has  now  been  quantized  and  is  in  a  suitable  form 
for  recursive  parameter  estimation  (Ref  2:319). 

The  joint  processing  algorithms  of  Eqs.  (6-23)  and 
(6-47)  were  derived  using  the  assumption  of  independent 
detector  noise  processes.  Development  of  an  algorithm  to 
include  the  case  of  spatially  correlated  noise  should  also 
be  considered.  Finally,  the  algorithm  of  Eq.  (6-47)  needs 
to  be  worked  out  to  an  implementable  form,  and  tested 
through  a  performance  evolution  by  means  of  Monte  Carlo 
simulation. 
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